##Install Packages if Needed
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("effectsize")) install.packages("effectsize")
if (!require("emmeans")) install.packages("emmeans")
if (!require("dplyr")) install.packages("dplyr")
if (!require("tidyr")) install.packages("tidyr")
if (!require("Rmisc")) install.packages("Rmisc")
if (!require("ggpubr")) install.packages("ggpubr")
if (!require("cowplot")) install.packages("cowplot")
if (!require("gridGraphics")) install.packages("gridGraphics")
if (!require("tidyr")) install.packages("tidyr")
##Load Packages
library(ggplot2) #Required for plotting
library(effectsize) #Required for eta_squared effect sizes
library(emmeans) #Required for pairwise comparisons
library(dplyr) #Required to seperate columns in dataframe
library(tidyr) #Required for data organization
library(Rmisc) #Required for summarySE for summary statistics
library(ggpubr) #Required for adding pairwise p-values to plots with stat_pvalue_manual
library(cowplot) #Required for arranging ggplots
library(gridGraphics) #Required for adding labels to arranged plots
library(tidyr) #Required for reshaping datafrom from a wide to long format.
Note: Run “Graphing Parameters” section from 01_ExperimentalSetup.R file
Note: Full Data with Bleaching Metrics and Color Scores created in 04_Models.R file
#Load Data
FullData<-read.csv("Outputs/FullData.csv", header=TRUE)
#Set factor variables
FullData$TimeP<-factor(FullData$TimeP, levels=c("W2", "M1", "M4"), ordered=TRUE)
FullData$Site<-factor(FullData$Site, levels=c("KL", "SS"), ordered=TRUE)
FullData$Genotype<-factor(FullData$Genotype, levels=c("AC10", "AC12", "AC8"), ordered=TRUE)
FullData$Treatment<-factor(FullData$Treatment, levels=c("Control", "Heat"), ordered=TRUE)
FullData$Treat<-factor(FullData$Treat, levels=c("C", "H"), ordered=TRUE)
FullData$Seas<-factor(FullData$Seas, levels=c("Summer1", "Summer2", "Winter"), ordered=TRUE)
##Control
FullData_C<-subset(FullData, Treat=="C")
##Heated
FullData_H<-subset(FullData, Treat=="H")
Calculating retention as proportion remaining relative to corresponding control levels (0-1) for each site and genotype at each timepoint.
##Calculate averages for Control Treatment for each Site, Genotype, Timepoint, and Season
names(FullData_C)
[1] "ID" "TimeP" "Site" "Genotype" "Treat" "Treatment"
[7] "Seas" "Set" "Score_Full" "Score_TP" "RandN" "SA_cm2"
[13] "Chl_ug.cm2" "Sym10.6_cm2"
FullData_C.a<-aggregate(FullData_C[,c(9:10, 13:14)], list(FullData_C$Site, FullData_C$Genotype, FullData_C$TimeP, FullData_C$Seas), mean, na.action = na.omit)
names(FullData_C.a)[1:4]<-c("Site", "Genotype", "TimeP", "Seas")
names(FullData_C.a)[5:8]<-paste(names(FullData_C)[c(9:10, 13:14)], "C", sep="_")
##Merge Control Averages with Heated Samples
#Merges by Site, Genotype, Timepoint, and Season
TolData<-merge(FullData_H, FullData_C.a, all.x=TRUE )
##Calculate Proportion Retained for each Bleaching Metric relative to Control
TolData$Score_Full.prop<-round((TolData$Score_Full/TolData$Score_Full_C), 4)
TolData$Score_TP.prop<-round((TolData$Score_TP/TolData$Score_TP_C), 4)
TolData$Chl.prop<-round((TolData$Chl_ug.cm2/TolData$Chl_ug.cm2_C), 4)
TolData$Sym.prop<-round((TolData$Sym10.6_cm2/TolData$Sym10.6_cm2_C), 4)
##Set values >1 to 1
TolData$Score_Full.prop[which(TolData$Score_Full.prop>1)]<-1.0000
TolData$Score_TP.prop[which(TolData$Score_TP.prop>1)]<-1.0000
TolData$Chl.prop[which(TolData$Chl.prop>1)]<-1.0000
TolData$Sym.prop[which(TolData$Sym.prop>1)]<-1.0000
##Create Site and Genotype Variable
TolData$Site.Geno<-paste(TolData$Site, TolData$Genotype, sep="_")
TolData_W2<-subset(TolData, TimeP=="W2")
TolData_M1<-subset(TolData, TimeP=="M1")
TolData_M4<-subset(TolData, TimeP=="M4")
##Check normality
hist(TolData$Chl.prop)
shapiro.test(TolData$Chl.prop)
Shapiro-Wilk normality test
data: TolData$Chl.prop
W = 0.97395, p-value = 0.165
#Normal
##Model as a function of Site and Genotype and Timepoint
Tol_Chl_lm<-lm(Chl.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Chl_lm)))
#Q-Q plot
qqnorm(resid(Tol_Chl_lm)); qqline(resid(Tol_Chl_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_lm), resid(Tol_Chl_lm))
#Model Results
summary(Tol_Chl_lm)
Call:
lm(formula = Chl.prop ~ Site + Genotype + TimeP + Site:Genotype +
Site:TimeP + Genotype:TimeP, data = TolData)
Residuals:
Min 1Q Median 3Q Max
-0.127060 -0.031006 -0.004219 0.036829 0.097690
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.2100753 0.0065828 31.913 < 2e-16 ***
Site.L -0.0093229 0.0092942 -1.003 0.320289
Genotype.L -0.0980285 0.0111500 -8.792 5.31e-12 ***
Genotype.Q -0.0467859 0.0116301 -4.023 0.000180 ***
TimeP.L 0.0638187 0.0112994 5.648 6.22e-07 ***
TimeP.Q 0.0464233 0.0115033 4.036 0.000173 ***
Site.L:Genotype.L 0.0003917 0.0157685 0.025 0.980272
Site.L:Genotype.Q 0.0592529 0.0164183 3.609 0.000673 ***
Site.L:TimeP.L -0.0528050 0.0159798 -3.304 0.001694 **
Site.L:TimeP.Q -0.0236749 0.0161894 -1.462 0.149435
Genotype.L:TimeP.L 0.0406686 0.0194339 2.093 0.041093 *
Genotype.Q:TimeP.L -0.0627603 0.0196595 -3.192 0.002354 **
Genotype.L:TimeP.Q 0.0800244 0.0191901 4.170 0.000111 ***
Genotype.Q:TimeP.Q -0.1328570 0.0206170 -6.444 3.28e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05393 on 54 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.8229, Adjusted R-squared: 0.7803
F-statistic: 19.3 on 13 and 54 DF, p-value: 8.612e-16
anova(Tol_Chl_lm)
Analysis of Variance Table
Response: Chl.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.001563 0.001563 0.5375 0.466633
Genotype 2 0.297307 0.148654 51.1117 3.495e-13 ***
TimeP 2 0.132285 0.066142 22.7418 6.845e-08 ***
Site:Genotype 2 0.041114 0.020557 7.0681 0.001877 **
Site:TimeP 2 0.038643 0.019322 6.6434 0.002634 **
Genotype:TimeP 4 0.218951 0.054738 18.8205 9.691e-10 ***
Residuals 54 0.157054 0.002908
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Chl_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
----------------------------------------
Site | 1.76e-03 | [0.00, 1.00]
Genotype | 0.34 | [0.16, 1.00]
TimeP | 0.15 | [0.02, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]
Site:TimeP | 0.04 | [0.00, 1.00]
Genotype:TimeP | 0.25 | [0.06, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_lm.res<-data.frame(anova(Tol_Chl_lm))
Tol_Chl_lm.res$Predictor<-rownames(Tol_Chl_lm.res)
Tol_Chl_lm.res$EtaSq<-c(eta_squared(Tol_Chl_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_lm.res))
Tol_Chl_lm.res<-Tol_Chl_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Strong influence of Timepoint, including interactions with main variables of interest, Site and Genotype. Will analyze each Timepoint individually.
##Check normality
hist(TolData_W2$Chl.prop)
shapiro.test(TolData_W2$Chl.prop)
Shapiro-Wilk normality test
data: TolData_W2$Chl.prop
W = 0.98407, p-value = 0.9671
#Normal
##Model as a function of Site and Genotype
Tol_Chl_W2_lm<-lm(Chl.prop~Site+Genotype+Site:Genotype, data=TolData_W2)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Chl_W2_lm)))
#Q-Q plot
qqnorm(resid(Tol_Chl_W2_lm)); qqline(resid(Tol_Chl_W2_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_W2_lm), resid(Tol_Chl_W2_lm))
#Model Results
summary(Tol_Chl_W2_lm)
Call:
lm(formula = Chl.prop ~ Site + Genotype + Site:Genotype, data = TolData_W2)
Residuals:
Min 1Q Median 3Q Max
-0.089900 -0.031587 -0.008183 0.038881 0.103550
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.1842264 0.0119450 15.423 5.03e-11 ***
Site.L 0.0188110 0.0168928 1.114 0.281921
Genotype.L -0.0947759 0.0204291 -4.639 0.000273 ***
Genotype.Q -0.0585870 0.0209464 -2.797 0.012921 *
Site.L:Genotype.L -0.0002833 0.0288911 -0.010 0.992297
Site.L:Genotype.Q 0.0366569 0.0296227 1.237 0.233769
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05552 on 16 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.6758, Adjusted R-squared: 0.5745
F-statistic: 6.671 on 5 and 16 DF, p-value: 0.001551
anova(Tol_Chl_W2_lm)
Analysis of Variance Table
Response: Chl.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.005580 0.005580 1.8106 0.197195
Genotype 2 0.092482 0.046241 15.0038 0.000214 ***
Site:Genotype 2 0.004732 0.002366 0.7677 0.480448
Residuals 16 0.049311 0.003082
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Chl_W2_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.04 | [0.00, 1.00]
Genotype | 0.61 | [0.29, 1.00]
Site:Genotype | 0.03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_W2_lm.res<-data.frame(anova(Tol_Chl_W2_lm))
Tol_Chl_W2_lm.res$Predictor<-rownames(Tol_Chl_W2_lm.res)
Tol_Chl_W2_lm.res$EtaSq<-c(eta_squared(Tol_Chl_W2_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_W2_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_W2_lm.res))
Tol_Chl_W2_lm.res$TimeP<-rep("W2", nrow(Tol_Chl_W2_lm.res))
Tol_Chl_W2_lm.res<-Tol_Chl_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Chl_W2_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.2033 0.0278 16 0.1445 0.262
AC12 0.2399 0.0278 16 0.1811 0.299
AC8 0.0696 0.0278 16 0.0107 0.128
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.2514 0.0278 16 0.1925 0.310
AC12 0.2242 0.0321 16 0.1563 0.292
AC8 0.1170 0.0321 16 0.0491 0.185
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 -0.0366 0.0393 16 -0.933 0.6280
AC10 - AC8 0.1338 0.0393 16 3.407 0.0095
AC12 - AC8 0.1704 0.0393 16 4.340 0.0014
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0272 0.0424 16 0.640 0.8003
AC10 - AC8 0.1343 0.0424 16 3.168 0.0156
AC12 - AC8 0.1072 0.0453 16 2.364 0.0753
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Chl_W2_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.2033 0.0278 16 0.1445 0.262
SS 0.2514 0.0278 16 0.1925 0.310
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.2399 0.0278 16 0.1811 0.299
SS 0.2242 0.0321 16 0.1563 0.292
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.0696 0.0278 16 0.0107 0.128
SS 0.1170 0.0321 16 0.0491 0.185
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0481 0.0393 16 -1.224 0.2387
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.0157 0.0424 16 0.371 0.7156
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.0475 0.0424 16 -1.120 0.2793
##Save p-values
#Genotypes within Sites
Tol_Chl_W2_lm.geno<-data.frame(emmeans(Tol_Chl_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_Chl_W2_lm.geno<-Tol_Chl_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_W2_lm.geno$group1<-paste(Tol_Chl_W2_lm.geno$Site, Tol_Chl_W2_lm.geno$group1, sep="_")
Tol_Chl_W2_lm.geno$group2<-paste(Tol_Chl_W2_lm.geno$Site, Tol_Chl_W2_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Chl_W2_lm.site<-data.frame(emmeans(Tol_Chl_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_Chl_W2_lm.site<-Tol_Chl_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_W2_lm.site$group1<-paste(Tol_Chl_W2_lm.site$group1, Tol_Chl_W2_lm.site$Genotype, sep="_")
Tol_Chl_W2_lm.site$group2<-paste(Tol_Chl_W2_lm.site$group2, Tol_Chl_W2_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Chl_W2_lm.p<-rbind(Tol_Chl_W2_lm.geno[,c(1:2,4:8)], Tol_Chl_W2_lm.site[,c(1:2,4:8)])
Tol_Chl_W2_lm.p<-Tol_Chl_W2_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Chl_W2_lm.p$Sig<-ifelse(Tol_Chl_W2_lm.p$p<0.001, "***", ifelse(Tol_Chl_W2_lm.p$p<0.01, "**", ifelse(Tol_Chl_W2_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Chl_W2_lm.p$Response<-rep("Chlorophyll", nrow(Tol_Chl_W2_lm.p))
Tol_Chl_W2_lm.p$TimeP<-rep("W2", nrow(Tol_Chl_W2_lm.p))
##Summary statistics by Site and Genotype
Tol_Chl_W2_SG<-summarySE(TolData_W2, measurevar="Chl.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Chl_W2_SG.plot<-ggplot(Tol_Chl_W2_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Chlorophyll Retained")+
ylim(0, 1)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Chl_W2_lm.p, y.position=0.4, step.increase=0.25, label="Sig", hide.ns=TRUE); Tol_Chl_W2_SG.plot
##Check normality
hist(TolData_M1$Chl.prop)
shapiro.test(TolData_M1$Chl.prop)
Shapiro-Wilk normality test
data: TolData_M1$Chl.prop
W = 0.9306, p-value = 0.1263
#Normal
##Model as a function of Site and Genotype
Tol_Chl_M1_lm<-lm(Chl.prop~Site+Genotype+Site:Genotype, data=TolData_M1)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Chl_M1_lm)))
#Q-Q plot
qqnorm(resid(Tol_Chl_M1_lm)); qqline(resid(Tol_Chl_M1_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_M1_lm), resid(Tol_Chl_M1_lm))
#Model Results
summary(Tol_Chl_M1_lm)
Call:
lm(formula = Chl.prop ~ Site + Genotype + Site:Genotype, data = TolData_M1)
Residuals:
Min 1Q Median 3Q Max
-0.107050 -0.034575 0.003479 0.025531 0.088050
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.17217 0.01178 14.619 1.12e-10 ***
Site.L 0.01036 0.01666 0.622 0.5428
Genotype.L -0.16337 0.01935 -8.442 2.74e-07 ***
Genotype.Q 0.06169 0.02139 2.883 0.0108 *
Site.L:Genotype.L 0.03819 0.02737 1.395 0.1820
Site.L:Genotype.Q 0.05454 0.03026 1.803 0.0903 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05474 on 16 degrees of freedom
Multiple R-squared: 0.8424, Adjusted R-squared: 0.7932
F-statistic: 17.11 on 5 and 16 DF, p-value: 6.411e-06
anova(Tol_Chl_M1_lm)
Analysis of Variance Table
Response: Chl.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.002283 0.002283 0.7619 0.3956
Genotype 2 0.238424 0.119212 39.7888 6.168e-07 ***
Site:Genotype 2 0.015569 0.007785 2.5982 0.1054
Residuals 16 0.047938 0.002996
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Chl_M1_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
---------------------------------------
Site | 7.50e-03 | [0.00, 1.00]
Genotype | 0.78 | [0.58, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_M1_lm.res<-data.frame(anova(Tol_Chl_M1_lm))
Tol_Chl_M1_lm.res$Predictor<-rownames(Tol_Chl_M1_lm.res)
Tol_Chl_M1_lm.res$EtaSq<-c(eta_squared(Tol_Chl_M1_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_M1_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_M1_lm.res))
Tol_Chl_M1_lm.res$TimeP<-rep("M1", nrow(Tol_Chl_M1_lm.res))
Tol_Chl_M1_lm.res<-Tol_Chl_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Chl_M1_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.3089 0.0274 16 0.2509 0.3669
AC12 0.1460 0.0316 16 0.0790 0.2130
AC8 0.0397 0.0274 16 -0.0183 0.0977
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.3169 0.0274 16 0.2588 0.3749
AC12 0.0976 0.0316 16 0.0306 0.1646
AC8 0.1240 0.0274 16 0.0660 0.1820
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.1629 0.0418 16 3.897 0.0035
AC10 - AC8 0.2692 0.0387 16 6.956 <.0001
AC12 - AC8 0.1063 0.0418 16 2.543 0.0538
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.2192 0.0418 16 5.244 0.0002
AC10 - AC8 0.1928 0.0387 16 4.983 0.0004
AC12 - AC8 -0.0264 0.0418 16 -0.631 0.8056
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Chl_M1_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.3089 0.0274 16 0.2509 0.3669
SS 0.3169 0.0274 16 0.2588 0.3749
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.1460 0.0316 16 0.0790 0.2130
SS 0.0976 0.0316 16 0.0306 0.1646
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.0397 0.0274 16 -0.0183 0.0977
SS 0.1240 0.0274 16 0.0660 0.1820
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.00795 0.0387 16 -0.205 0.8398
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.04833 0.0447 16 1.081 0.2955
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.08432 0.0387 16 -2.179 0.0447
##Save p-values
#Genotypes within Sites
Tol_Chl_M1_lm.geno<-data.frame(emmeans(Tol_Chl_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_Chl_M1_lm.geno<-Tol_Chl_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M1_lm.geno$group1<-paste(Tol_Chl_M1_lm.geno$Site, Tol_Chl_M1_lm.geno$group1, sep="_")
Tol_Chl_M1_lm.geno$group2<-paste(Tol_Chl_M1_lm.geno$Site, Tol_Chl_M1_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Chl_M1_lm.site<-data.frame(emmeans(Tol_Chl_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_Chl_M1_lm.site<-Tol_Chl_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M1_lm.site$group1<-paste(Tol_Chl_M1_lm.site$group1, Tol_Chl_M1_lm.site$Genotype, sep="_")
Tol_Chl_M1_lm.site$group2<-paste(Tol_Chl_M1_lm.site$group2, Tol_Chl_M1_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Chl_M1_lm.p<-rbind(Tol_Chl_M1_lm.geno[,c(1:2,4:8)], Tol_Chl_M1_lm.site[,c(1:2,4:8)])
Tol_Chl_M1_lm.p<-Tol_Chl_M1_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Chl_M1_lm.p$Sig<-ifelse(Tol_Chl_M1_lm.p$p<0.001, "***", ifelse(Tol_Chl_M1_lm.p$p<0.01, "**", ifelse(Tol_Chl_M1_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Chl_M1_lm.p$Response<-rep("Chlorophyll", nrow(Tol_Chl_M1_lm.p))
Tol_Chl_M1_lm.p$TimeP<-rep("M1", nrow(Tol_Chl_M1_lm.p))
##Summary statistics by Site and Genotype
Tol_Chl_M1_SG<-summarySE(TolData_M1, measurevar="Chl.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Chl_M1_SG.plot<-ggplot(Tol_Chl_M1_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Chlorophyll Retained")+
ylim(0, 1)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Chl_M1_lm.p, y.position=0.4, step.increase=0.25, label="Sig", hide.ns=TRUE); Tol_Chl_M1_SG.plot
##Check normality
hist(TolData_M4$Chl.prop)
shapiro.test(TolData_M4$Chl.prop)
Shapiro-Wilk normality test
data: TolData_M4$Chl.prop
W = 0.91299, p-value = 0.04096
#Not Normal
##Try log+1 transformation
hist(log(TolData_M4$Chl.prop+1))
shapiro.test(log(TolData_M4$Chl.prop+1))
Shapiro-Wilk normality test
data: log(TolData_M4$Chl.prop + 1)
W = 0.93102, p-value = 0.1028
#Normal
##Model as a function of Site and Genotype
##Model with log transformation
Tol_Chl_M4_lm<-lm(log(Chl.prop+1)~Site+Genotype+Site:Genotype, data=TolData_M4)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Chl_M4_lm)))
#Q-Q plot
qqnorm(resid(Tol_Chl_M4_lm)); qqline(resid(Tol_Chl_M4_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_M4_lm), resid(Tol_Chl_M4_lm))
#Model Results
summary(Tol_Chl_M4_lm)
Call:
lm(formula = log(Chl.prop + 1) ~ Site + Genotype + Site:Genotype,
data = TolData_M4)
Residuals:
Min 1Q Median 3Q Max
-0.062753 -0.025900 -0.004896 0.017009 0.065759
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.23863 0.00807 29.568 < 2e-16 ***
Site.L -0.04272 0.01141 -3.743 0.00149 **
Genotype.L -0.03027 0.01398 -2.166 0.04400 *
Genotype.Q -0.11052 0.01398 -7.907 2.9e-07 ***
Site.L:Genotype.L -0.03189 0.01977 -1.613 0.12415
Site.L:Genotype.Q 0.05753 0.01977 2.910 0.00934 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.03954 on 18 degrees of freedom
Multiple R-squared: 0.8368, Adjusted R-squared: 0.7915
F-statistic: 18.46 on 5 and 18 DF, p-value: 1.598e-06
anova(Tol_Chl_M4_lm)
Analysis of Variance Table
Response: log(Chl.prop + 1)
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.021901 0.021901 14.0102 0.001489 **
Genotype 2 0.105056 0.052528 33.6031 8.379e-07 ***
Site:Genotype 2 0.017304 0.008652 5.5349 0.013381 *
Residuals 18 0.028137 0.001563
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Chl_M4_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.13 | [0.00, 1.00]
Genotype | 0.61 | [0.32, 1.00]
Site:Genotype | 0.10 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_M4_lm.res<-data.frame(anova(Tol_Chl_M4_lm))
Tol_Chl_M4_lm.res$Predictor<-rownames(Tol_Chl_M4_lm.res)
Tol_Chl_M4_lm.res$EtaSq<-c(eta_squared(Tol_Chl_M4_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_M4_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_M4_lm.res))
Tol_Chl_M4_lm.res$TimeP<-rep("M4", nrow(Tol_Chl_M4_lm.res))
Tol_Chl_M4_lm.res<-Tol_Chl_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype and Site*Genotype.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Chl_M4_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.213 0.0198 18 0.171 0.254
AC12 0.392 0.0198 18 0.351 0.434
AC8 0.202 0.0198 18 0.160 0.243
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.217 0.0198 18 0.176 0.259
AC12 0.265 0.0198 18 0.224 0.307
AC8 0.143 0.0198 18 0.101 0.184
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 -0.1797 0.028 18 -6.428 <.0001
AC10 - AC8 0.0109 0.028 18 0.391 0.9196
AC12 - AC8 0.1906 0.028 18 6.819 <.0001
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 -0.0482 0.028 18 -1.724 0.2237
AC10 - AC8 0.0747 0.028 18 2.672 0.0393
AC12 - AC8 0.1229 0.028 18 4.396 0.0010
Note: contrasts are still on the log(mu + 1) scale
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Chl_M4_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.213 0.0198 18 0.171 0.254
SS 0.217 0.0198 18 0.176 0.259
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.392 0.0198 18 0.351 0.434
SS 0.265 0.0198 18 0.224 0.307
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.202 0.0198 18 0.160 0.243
SS 0.143 0.0198 18 0.101 0.184
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.00468 0.028 18 -0.167 0.8688
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.12684 0.028 18 4.537 0.0003
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS 0.05909 0.028 18 2.114 0.0488
Note: contrasts are still on the log(mu + 1) scale
##Save p-values
#Genotypes within Sites
Tol_Chl_M4_lm.geno<-data.frame(emmeans(Tol_Chl_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_Chl_M4_lm.geno<-Tol_Chl_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M4_lm.geno$group1<-paste(Tol_Chl_M4_lm.geno$Site, Tol_Chl_M4_lm.geno$group1, sep="_")
Tol_Chl_M4_lm.geno$group2<-paste(Tol_Chl_M4_lm.geno$Site, Tol_Chl_M4_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Chl_M4_lm.site<-data.frame(emmeans(Tol_Chl_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_Chl_M4_lm.site<-Tol_Chl_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M4_lm.site$group1<-paste(Tol_Chl_M4_lm.site$group1, Tol_Chl_M4_lm.site$Genotype, sep="_")
Tol_Chl_M4_lm.site$group2<-paste(Tol_Chl_M4_lm.site$group2, Tol_Chl_M4_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Chl_M4_lm.p<-rbind(Tol_Chl_M4_lm.geno[,c(1:2,4:8)], Tol_Chl_M4_lm.site[,c(1:2,4:8)])
Tol_Chl_M4_lm.p<-Tol_Chl_M4_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Chl_M4_lm.p$Sig<-ifelse(Tol_Chl_M4_lm.p$p<0.001, "***", ifelse(Tol_Chl_M4_lm.p$p<0.01, "**", ifelse(Tol_Chl_M4_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Chl_M4_lm.p$Response<-rep("Chlorophyll", nrow(Tol_Chl_M4_lm.p))
Tol_Chl_M4_lm.p$TimeP<-rep("M4", nrow(Tol_Chl_M4_lm.p))
##Summary statistics by Site and Genotype
Tol_Chl_M4_SG<-summarySE(TolData_M4, measurevar="Chl.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Chl_M4_SG.plot<-ggplot(Tol_Chl_M4_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Chlorophyll Retained")+
ylim(0, 1)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Chl_M4_lm.p, y.position=0.55, step.increase=0.20, label="Sig", hide.ns=TRUE); Tol_Chl_M4_SG.plot
##Check normality
hist(TolData$Sym.prop)
shapiro.test(TolData$Sym.prop)
Shapiro-Wilk normality test
data: TolData$Sym.prop
W = 0.92235, p-value = 0.0003707
#Not Normal
##Try log+1 transformation
hist(log(TolData$Sym.prop+1))
shapiro.test(log(TolData$Sym.prop+1))
Shapiro-Wilk normality test
data: log(TolData$Sym.prop + 1)
W = 0.92164, p-value = 0.000345
#Not normal
##Try square transformation
hist((TolData$Sym.prop)^2)
shapiro.test((TolData$Sym.prop)^2)
Shapiro-Wilk normality test
data: (TolData$Sym.prop)^2
W = 0.87147, p-value = 3.933e-06
#Not normal
##Try cubed transformation
hist((TolData$Sym.prop)^3)
shapiro.test((TolData$Sym.prop)^3)
Shapiro-Wilk normality test
data: (TolData$Sym.prop)^3
W = 0.80697, p-value = 4.342e-08
#Not normal
##Model as a function of Site and Genotype and Timepoint
##Model with no transformation and check residuals
Tol_Sym_lm<-lm(Sym.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Sym_lm)))
#Q-Q plot
qqnorm(resid(Tol_Sym_lm)); qqline(resid(Tol_Sym_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_lm), resid(Tol_Sym_lm))
#Model Results
summary(Tol_Sym_lm)
Call:
lm(formula = Sym.prop ~ Site + Genotype + TimeP + Site:Genotype +
Site:TimeP + Genotype:TimeP, data = TolData)
Residuals:
Min 1Q Median 3Q Max
-0.37464 -0.08630 -0.00192 0.08972 0.35883
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.544513 0.018576 29.313 < 2e-16 ***
Site.L -0.001834 0.026197 -0.070 0.944443
Genotype.L -0.373179 0.031737 -11.758 < 2e-16 ***
Genotype.Q 0.028101 0.032605 0.862 0.392497
TimeP.L 0.133503 0.031737 4.206 9.64e-05 ***
TimeP.Q 0.190344 0.032605 5.838 2.94e-07 ***
Site.L:Genotype.L -0.060389 0.044883 -1.345 0.183999
Site.L:Genotype.Q 0.074992 0.045930 1.633 0.108240
Site.L:TimeP.L -0.159707 0.044883 -3.558 0.000778 ***
Site.L:TimeP.Q -0.027928 0.045930 -0.608 0.545653
Genotype.L:TimeP.L 0.085371 0.055316 1.543 0.128489
Genotype.Q:TimeP.L -0.084924 0.054623 -1.555 0.125749
Genotype.L:TimeP.Q 0.065986 0.054623 1.208 0.232204
Genotype.Q:TimeP.Q -0.205029 0.058264 -3.519 0.000878 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1535 on 55 degrees of freedom
Multiple R-squared: 0.8036, Adjusted R-squared: 0.7572
F-statistic: 17.31 on 13 and 55 DF, p-value: 6.022e-15
anova(Tol_Sym_lm)
Analysis of Variance Table
Response: Sym.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.0007 0.00068 0.0289 0.865708
Genotype 2 3.2796 1.63981 69.5880 8.602e-16 ***
TimeP 2 1.1556 0.57780 24.5197 2.438e-08 ***
Site:Genotype 2 0.1111 0.05557 2.3581 0.104096
Site:TimeP 2 0.3138 0.15688 6.6576 0.002575 **
Genotype:TimeP 4 0.4428 0.11069 4.6974 0.002474 **
Residuals 55 1.2961 0.02356
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Sym_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
----------------------------------------
Site | 1.03e-04 | [0.00, 1.00]
Genotype | 0.50 | [0.33, 1.00]
TimeP | 0.18 | [0.04, 1.00]
Site:Genotype | 0.02 | [0.00, 1.00]
Site:TimeP | 0.05 | [0.00, 1.00]
Genotype:TimeP | 0.07 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_lm.res<-data.frame(anova(Tol_Sym_lm))
Tol_Sym_lm.res$Predictor<-rownames(Tol_Sym_lm.res)
Tol_Sym_lm.res$EtaSq<-c(eta_squared(Tol_Sym_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_lm.res))
Tol_Sym_lm.res<-Tol_Sym_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Strong influence of Timepoint, including interactions with main variables of interest, Site and Genotype. Will analyze each Timepoint individually.
##Check normality
hist(TolData_W2$Sym.prop)
shapiro.test(TolData_W2$Sym.prop)
Shapiro-Wilk normality test
data: TolData_W2$Sym.prop
W = 0.93605, p-value = 0.1477
#Normal
##Model as a function of Site and Genotype
Tol_Sym_W2_lm<-lm(Sym.prop~Site+Genotype+Site:Genotype, data=TolData_W2)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Sym_W2_lm)))
#Q-Q plot
qqnorm(resid(Tol_Sym_W2_lm)); qqline(resid(Tol_Sym_W2_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_W2_lm), resid(Tol_Sym_W2_lm))
#Model Results
summary(Tol_Sym_W2_lm)
Call:
lm(formula = Sym.prop ~ Site + Genotype + Site:Genotype, data = TolData_W2)
Residuals:
Min 1Q Median 3Q Max
-0.26060 -0.11145 0.00000 0.08525 0.28580
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5237625 0.0325072 16.112 9.90e-12 ***
Site.L 0.0939568 0.0459722 2.044 0.0568 .
Genotype.L -0.4152131 0.0570402 -7.279 1.29e-06 ***
Genotype.Q -0.0005205 0.0555584 -0.009 0.9926
Site.L:Genotype.L -0.2061000 0.0806671 -2.555 0.0205 *
Site.L:Genotype.Q 0.0462891 0.0785715 0.589 0.5635
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.155 on 17 degrees of freedom
Multiple R-squared: 0.7924, Adjusted R-squared: 0.7313
F-statistic: 12.98 on 5 and 17 DF, p-value: 2.647e-05
anova(Tol_Sym_W2_lm)
Analysis of Variance Table
Response: Sym.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.16148 0.16148 6.7208 0.01897 *
Genotype 2 1.22843 0.61421 25.5640 7.508e-06 ***
Site:Genotype 2 0.16883 0.08441 3.5133 0.05283 .
Residuals 17 0.40845 0.02403
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Sym_W2_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.08 | [0.00, 1.00]
Genotype | 0.62 | [0.33, 1.00]
Site:Genotype | 0.09 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_W2_lm.res<-data.frame(anova(Tol_Sym_W2_lm))
Tol_Sym_W2_lm.res$Predictor<-rownames(Tol_Sym_W2_lm.res)
Tol_Sym_W2_lm.res$EtaSq<-c(eta_squared(Tol_Sym_W2_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_W2_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_W2_lm.res))
Tol_Sym_W2_lm.res$TimeP<-rep("W2", nrow(Tol_Sym_W2_lm.res))
Tol_Sym_W2_lm.res<-Tol_Sym_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Marginal effect (p<0.1) of Site * Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Sym_W2_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.634 0.0775 17 0.4708 0.798
AC12 0.484 0.0775 17 0.3210 0.648
AC8 0.253 0.0775 17 0.0897 0.417
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 1.000 0.0775 17 0.8365 1.164
AC12 0.564 0.0775 17 0.4004 0.727
AC8 0.207 0.0895 17 0.0179 0.396
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.150 0.110 17 1.367 0.3795
AC10 - AC8 0.381 0.110 17 3.477 0.0077
AC12 - AC8 0.231 0.110 17 2.110 0.1174
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.436 0.110 17 3.979 0.0026
AC10 - AC8 0.793 0.118 17 6.701 <.0001
AC12 - AC8 0.357 0.118 17 3.017 0.0201
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Sym_W2_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.634 0.0775 17 0.4708 0.798
SS 1.000 0.0775 17 0.8365 1.164
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.484 0.0775 17 0.3210 0.648
SS 0.564 0.0775 17 0.4004 0.727
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.253 0.0775 17 0.0897 0.417
SS 0.207 0.0895 17 0.0179 0.396
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.3657 0.110 17 -3.337 0.0039
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.0794 0.110 17 -0.725 0.4785
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS 0.0465 0.118 17 0.393 0.6994
##Save p-values
#Genotypes within Sites
Tol_Sym_W2_lm.geno<-data.frame(emmeans(Tol_Sym_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_Sym_W2_lm.geno<-Tol_Sym_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_W2_lm.geno$group1<-paste(Tol_Sym_W2_lm.geno$Site, Tol_Sym_W2_lm.geno$group1, sep="_")
Tol_Sym_W2_lm.geno$group2<-paste(Tol_Sym_W2_lm.geno$Site, Tol_Sym_W2_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Sym_W2_lm.site<-data.frame(emmeans(Tol_Sym_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_Sym_W2_lm.site<-Tol_Sym_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_W2_lm.site$group1<-paste(Tol_Sym_W2_lm.site$group1, Tol_Sym_W2_lm.site$Genotype, sep="_")
Tol_Sym_W2_lm.site$group2<-paste(Tol_Sym_W2_lm.site$group2, Tol_Sym_W2_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Sym_W2_lm.p<-rbind(Tol_Sym_W2_lm.geno[,c(1:2,4:8)], Tol_Sym_W2_lm.site[,c(1:2,4:8)])
Tol_Sym_W2_lm.p<-Tol_Sym_W2_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Sym_W2_lm.p$Sig<-ifelse(Tol_Sym_W2_lm.p$p<0.001, "***", ifelse(Tol_Sym_W2_lm.p$p<0.01, "**", ifelse(Tol_Sym_W2_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Sym_W2_lm.p$Response<-rep("Symbionts", nrow(Tol_Sym_W2_lm.p))
Tol_Sym_W2_lm.p$TimeP<-rep("W2", nrow(Tol_Sym_W2_lm.p))
##Summary statistics by Site and Genotype
Tol_Sym_W2_SG<-summarySE(TolData_W2, measurevar="Sym.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Sym_W2_SG.plot<-ggplot(Tol_Sym_W2_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Symbionts Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Sym_W2_lm.p, y.position=0.95, step.increase=0.15, label="Sig", hide.ns=TRUE); Tol_Sym_W2_SG.plot
##Check normality
hist(TolData_M1$Sym.prop)
shapiro.test(TolData_M1$Sym.prop)
Shapiro-Wilk normality test
data: TolData_M1$Sym.prop
W = 0.85748, p-value = 0.004618
#Not normal
##Try log+1 transformation
hist(log(TolData_M1$Sym.prop+1))
shapiro.test(log(TolData_M1$Sym.prop+1))
Shapiro-Wilk normality test
data: log(TolData_M1$Sym.prop + 1)
W = 0.87335, p-value = 0.009041
#Not normal but improved
##Model as a function of Site and Genotype
##Model with log transformation
Tol_Sym_M1_lm<-lm(log(Sym.prop+1)~Site+Genotype+Site:Genotype, data=TolData_M1)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Sym_M1_lm)))
#Q-Q plot
qqnorm(resid(Tol_Sym_M1_lm)); qqline(resid(Tol_Sym_M1_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_M1_lm), resid(Tol_Sym_M1_lm))
#Model Results
summary(Tol_Sym_M1_lm)
Call:
lm(formula = log(Sym.prop + 1) ~ Site + Genotype + Site:Genotype,
data = TolData_M1)
Residuals:
Min 1Q Median 3Q Max
-0.098079 -0.055799 -0.004888 0.043331 0.151113
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.307716 0.016556 18.587 2.95e-12 ***
Site.L 0.027767 0.023413 1.186 0.252954
Genotype.L -0.295494 0.027204 -10.862 8.59e-09 ***
Genotype.Q 0.129387 0.030075 4.302 0.000548 ***
Site.L:Genotype.L 0.080682 0.038472 2.097 0.052224 .
Site.L:Genotype.Q -0.005572 0.042532 -0.131 0.897412
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.07694 on 16 degrees of freedom
Multiple R-squared: 0.8989, Adjusted R-squared: 0.8673
F-statistic: 28.46 on 5 and 16 DF, p-value: 2e-07
anova(Tol_Sym_M1_lm)
Analysis of Variance Table
Response: log(Sym.prop + 1)
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.00823 0.00823 1.3902 0.2556
Genotype 2 0.80811 0.40406 68.2486 1.468e-08 ***
Site:Genotype 2 0.02614 0.01307 2.2076 0.1423
Residuals 16 0.09473 0.00592
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Sym_M1_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
---------------------------------------
Site | 8.78e-03 | [0.00, 1.00]
Genotype | 0.86 | [0.73, 1.00]
Site:Genotype | 0.03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_M1_lm.res<-data.frame(anova(Tol_Sym_M1_lm))
Tol_Sym_M1_lm.res$Predictor<-rownames(Tol_Sym_M1_lm.res)
Tol_Sym_M1_lm.res$EtaSq<-c(eta_squared(Tol_Sym_M1_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_M1_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_M1_lm.res))
Tol_Sym_M1_lm.res$TimeP<-rep("M1", nrow(Tol_Sym_M1_lm.res))
Tol_Sym_M1_lm.res<-Tol_Sym_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Sym_M1_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.5918 0.0385 16 0.5102 0.673
AC12 0.1792 0.0444 16 0.0850 0.273
AC8 0.0932 0.0385 16 0.0117 0.175
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.5472 0.0385 16 0.4656 0.629
AC12 0.2249 0.0444 16 0.1307 0.319
AC8 0.2100 0.0385 16 0.1284 0.292
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.413 0.0588 16 7.021 <.0001
AC10 - AC8 0.499 0.0544 16 9.164 <.0001
AC12 - AC8 0.086 0.0588 16 1.463 0.3339
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.322 0.0588 16 5.483 0.0001
AC10 - AC8 0.337 0.0544 16 6.198 <.0001
AC12 - AC8 0.015 0.0588 16 0.255 0.9650
Note: contrasts are still on the log(mu + 1) scale
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Sym_M1_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.5918 0.0385 16 0.5102 0.673
SS 0.5472 0.0385 16 0.4656 0.629
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.1792 0.0444 16 0.0850 0.273
SS 0.2249 0.0444 16 0.1307 0.319
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.0932 0.0385 16 0.0117 0.175
SS 0.2100 0.0385 16 0.1284 0.292
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS 0.0446 0.0544 16 0.820 0.4241
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.0457 0.0628 16 -0.727 0.4775
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.1167 0.0544 16 -2.146 0.0476
Note: contrasts are still on the log(mu + 1) scale
##Save p-values
#Genotypes within Sites
Tol_Sym_M1_lm.geno<-data.frame(emmeans(Tol_Sym_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_Sym_M1_lm.geno<-Tol_Sym_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M1_lm.geno$group1<-paste(Tol_Sym_M1_lm.geno$Site, Tol_Sym_M1_lm.geno$group1, sep="_")
Tol_Sym_M1_lm.geno$group2<-paste(Tol_Sym_M1_lm.geno$Site, Tol_Sym_M1_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Sym_M1_lm.site<-data.frame(emmeans(Tol_Sym_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_Sym_M1_lm.site<-Tol_Sym_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M1_lm.site$group1<-paste(Tol_Sym_M1_lm.site$group1, Tol_Sym_M1_lm.site$Genotype, sep="_")
Tol_Sym_M1_lm.site$group2<-paste(Tol_Sym_M1_lm.site$group2, Tol_Sym_M1_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Sym_M1_lm.p<-rbind(Tol_Sym_M1_lm.geno[,c(1:2,4:8)], Tol_Sym_M1_lm.site[,c(1:2,4:8)])
Tol_Sym_M1_lm.p<-Tol_Sym_M1_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Sym_M1_lm.p$Sig<-ifelse(Tol_Sym_M1_lm.p$p<0.001, "***", ifelse(Tol_Sym_M1_lm.p$p<0.01, "**", ifelse(Tol_Sym_M1_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Sym_M1_lm.p$Response<-rep("Symbionts", nrow(Tol_Sym_M1_lm.p))
Tol_Sym_M1_lm.p$TimeP<-rep("M1", nrow(Tol_Sym_M1_lm.p))
##Summary statistics by Site and Genotype
Tol_Sym_M1_SG<-summarySE(TolData_M1, measurevar="Sym.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Sym_M1_SG.plot<-ggplot(Tol_Sym_M1_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Symbionts Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Sym_M1_lm.p, y.position=0.95, step.increase=0.15, label="Sig", hide.ns=TRUE); Tol_Sym_M1_SG.plot
##Check normality
hist(TolData_M4$Sym.prop)
shapiro.test(TolData_M4$Sym.prop)
Shapiro-Wilk normality test
data: TolData_M4$Sym.prop
W = 0.88557, p-value = 0.01078
#Not Normal
##Try square transformation
hist((TolData_M4$Sym.prop)^2)
shapiro.test((TolData_M4$Sym.prop)^2)
Shapiro-Wilk normality test
data: (TolData_M4$Sym.prop)^2
W = 0.87973, p-value = 0.008208
#Not Normal
##Try cubed transformation
hist((TolData_M4$Sym.prop)^3)
shapiro.test((TolData_M4$Sym.prop)^3)
Shapiro-Wilk normality test
data: (TolData_M4$Sym.prop)^3
W = 0.8488, p-value = 0.002077
#Not Normal
##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_Sym_M4_lm<-lm(Sym.prop~Site+Genotype+Site:Genotype, data=TolData_M4)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Sym_M4_lm)))
#Q-Q plot
qqnorm(resid(Tol_Sym_M4_lm)); qqline(resid(Tol_Sym_M4_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_M4_lm), resid(Tol_Sym_M4_lm))
#Model Results
summary(Tol_Sym_M4_lm)
Call:
lm(formula = Sym.prop ~ Site + Genotype + Site:Genotype, data = TolData_M4)
Residuals:
Min 1Q Median 3Q Max
-0.31820 -0.08804 0.00000 0.08618 0.29740
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.71662 0.03157 22.702 1.07e-14 ***
Site.L -0.12617 0.04464 -2.826 0.0112 *
Genotype.L -0.28587 0.05467 -5.229 5.67e-05 ***
Genotype.Q -0.11565 0.05467 -2.115 0.0486 *
Site.L:Genotype.L -0.09494 0.07732 -1.228 0.2353
Site.L:Genotype.Q 0.17275 0.07732 2.234 0.0384 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1546 on 18 degrees of freedom
Multiple R-squared: 0.7201, Adjusted R-squared: 0.6423
F-statistic: 9.26 on 5 and 18 DF, p-value: 0.0001686
anova(Tol_Sym_M4_lm)
Analysis of Variance Table
Response: Sym.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.19101 0.19101 7.9874 0.011190 *
Genotype 2 0.76080 0.38040 15.9068 0.000105 ***
Site:Genotype 2 0.15542 0.07771 3.2496 0.062386 .
Residuals 18 0.43046 0.02391
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Sym_M4_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.12 | [0.00, 1.00]
Genotype | 0.49 | [0.18, 1.00]
Site:Genotype | 0.10 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_M4_lm.res<-data.frame(anova(Tol_Sym_M4_lm))
Tol_Sym_M4_lm.res$Predictor<-rownames(Tol_Sym_M4_lm.res)
Tol_Sym_M4_lm.res$EtaSq<-c(eta_squared(Tol_Sym_M4_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_M4_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_M4_lm.res))
Tol_Sym_M4_lm.res$TimeP<-rep("M4", nrow(Tol_Sym_M4_lm.res))
Tol_Sym_M4_lm.res<-Tol_Sym_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Marginal effect (p<0.1) of Site * Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Sym_M4_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.863 0.0773 18 0.701 1.026
AC12 1.000 0.0773 18 0.838 1.162
AC8 0.554 0.0773 18 0.392 0.717
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.880 0.0773 18 0.717 1.042
AC12 0.622 0.0773 18 0.460 0.785
AC8 0.380 0.0773 18 0.218 0.543
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 -0.137 0.109 18 -1.249 0.4409
AC10 - AC8 0.309 0.109 18 2.829 0.0285
AC12 - AC8 0.446 0.109 18 4.078 0.0019
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.258 0.109 18 2.356 0.0734
AC10 - AC8 0.499 0.109 18 4.565 0.0007
AC12 - AC8 0.242 0.109 18 2.210 0.0965
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Sym_M4_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.863 0.0773 18 0.701 1.026
SS 0.880 0.0773 18 0.717 1.042
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 1.000 0.0773 18 0.838 1.162
SS 0.622 0.0773 18 0.460 0.785
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.554 0.0773 18 0.392 0.717
SS 0.380 0.0773 18 0.218 0.543
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0163 0.109 18 -0.149 0.8835
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.3779 0.109 18 3.456 0.0028
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS 0.1736 0.109 18 1.588 0.1297
##Save p-values
#Genotypes within Sites
Tol_Sym_M4_lm.geno<-data.frame(emmeans(Tol_Sym_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_Sym_M4_lm.geno<-Tol_Sym_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M4_lm.geno$group1<-paste(Tol_Sym_M4_lm.geno$Site, Tol_Sym_M4_lm.geno$group1, sep="_")
Tol_Sym_M4_lm.geno$group2<-paste(Tol_Sym_M4_lm.geno$Site, Tol_Sym_M4_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Sym_M4_lm.site<-data.frame(emmeans(Tol_Sym_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_Sym_M4_lm.site<-Tol_Sym_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M4_lm.site$group1<-paste(Tol_Sym_M4_lm.site$group1, Tol_Sym_M4_lm.site$Genotype, sep="_")
Tol_Sym_M4_lm.site$group2<-paste(Tol_Sym_M4_lm.site$group2, Tol_Sym_M4_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Sym_M4_lm.p<-rbind(Tol_Sym_M4_lm.geno[,c(1:2,4:8)], Tol_Sym_M4_lm.site[,c(1:2,4:8)])
Tol_Sym_M4_lm.p<-Tol_Sym_M4_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Sym_M4_lm.p$Sig<-ifelse(Tol_Sym_M4_lm.p$p<0.001, "***", ifelse(Tol_Sym_M4_lm.p$p<0.01, "**", ifelse(Tol_Sym_M4_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Sym_M4_lm.p$Response<-rep("Symbionts", nrow(Tol_Sym_M4_lm.p))
Tol_Sym_M4_lm.p$TimeP<-rep("M4", nrow(Tol_Sym_M4_lm.p))
##Summary statistics by Site and Genotype
Tol_Sym_M4_SG<-summarySE(TolData_M4, measurevar="Sym.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Sym_M4_SG.plot<-ggplot(Tol_Sym_M4_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Symbionts Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Sym_M4_lm.p, y.position=1.1, step.increase=0.15, label="Sig", hide.ns=TRUE); Tol_Sym_M4_SG.plot
##Check normality
hist(TolData$Score_Full.prop)
shapiro.test(TolData$Score_Full.prop)
Shapiro-Wilk normality test
data: TolData$Score_Full.prop
W = 0.93829, p-value = 0.002023
#Not Normal
##Try square transformation
hist((TolData$Score_Full.prop)^2)
shapiro.test((TolData$Score_Full.prop)^2)
Shapiro-Wilk normality test
data: (TolData$Score_Full.prop)^2
W = 0.95035, p-value = 0.008158
#Not normal
##Try cubed transformation
hist((TolData$Score_Full.prop)^3)
shapiro.test((TolData$Score_Full.prop)^3)
Shapiro-Wilk normality test
data: (TolData$Score_Full.prop)^3
W = 0.95927, p-value = 0.02433
#Not normal
##Model as a function of Site and Genotype and Timepoint
##Model with no transformation and check residuals
Tol_ScoreF_lm<-lm(Score_Full.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreF_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreF_lm)); qqline(resid(Tol_ScoreF_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_lm), resid(Tol_ScoreF_lm))
#Model Results
summary(Tol_ScoreF_lm)
Call:
lm(formula = Score_Full.prop ~ Site + Genotype + TimeP + Site:Genotype +
Site:TimeP + Genotype:TimeP, data = TolData)
Residuals:
Min 1Q Median 3Q Max
-0.110993 -0.016292 0.006243 0.019125 0.105407
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.915091 0.004407 207.631 < 2e-16 ***
Site.L 0.010360 0.006215 1.667 0.101242
Genotype.L -0.066874 0.007530 -8.881 3.28e-12 ***
Genotype.Q 0.012290 0.007736 1.589 0.117866
TimeP.L 0.030163 0.007530 4.006 0.000187 ***
TimeP.Q 0.035927 0.007736 4.644 2.17e-05 ***
Site.L:Genotype.L 0.008168 0.010649 0.767 0.446384
Site.L:Genotype.Q 0.016902 0.010898 1.551 0.126651
Site.L:TimeP.L -0.014751 0.010649 -1.385 0.171595
Site.L:TimeP.Q -0.004576 0.010898 -0.420 0.676193
Genotype.L:TimeP.L 0.014217 0.013125 1.083 0.283412
Genotype.Q:TimeP.L 0.015151 0.012960 1.169 0.247420
Genotype.L:TimeP.Q 0.006442 0.012960 0.497 0.621134
Genotype.Q:TimeP.Q -0.048959 0.013824 -3.542 0.000819 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.03642 on 55 degrees of freedom
Multiple R-squared: 0.717, Adjusted R-squared: 0.6501
F-statistic: 10.72 on 13 and 55 DF, p-value: 7.67e-11
anova(Tol_ScoreF_lm)
Analysis of Variance Table
Response: Score_Full.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.004778 0.004778 3.6019 0.062967 .
Genotype 2 0.105770 0.052885 39.8672 1.992e-11 ***
TimeP 2 0.046815 0.023407 17.6455 1.202e-06 ***
Site:Genotype 2 0.004264 0.002132 1.6071 0.209738
Site:TimeP 2 0.003044 0.001522 1.1474 0.324936
Genotype:TimeP 4 0.020161 0.005040 3.7996 0.008476 **
Residuals 55 0.072959 0.001327
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreF_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
------------------------------------
Site | 0.02 | [0.00, 1.00]
Genotype | 0.41 | [0.24, 1.00]
TimeP | 0.18 | [0.04, 1.00]
Site:Genotype | 0.02 | [0.00, 1.00]
Site:TimeP | 0.01 | [0.00, 1.00]
Genotype:TimeP | 0.08 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreF_lm.res<-data.frame(anova(Tol_ScoreF_lm))
Tol_ScoreF_lm.res$Predictor<-rownames(Tol_ScoreF_lm.res)
Tol_ScoreF_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_lm.res))
Tol_ScoreF_lm.res<-Tol_ScoreF_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Strong influence of Timepoint, including interactions with Genotype. Will analyze each Timepoint individually.
##Check normality
hist(TolData_W2$Score_Full.prop)
shapiro.test(TolData_W2$Score_Full.prop)
Shapiro-Wilk normality test
data: TolData_W2$Score_Full.prop
W = 0.9067, p-value = 0.03483
#Not normal
##Try square transformation
hist((TolData_W2$Score_Full.prop)^2)
shapiro.test((TolData_W2$Score_Full.prop)^2)
Shapiro-Wilk normality test
data: (TolData_W2$Score_Full.prop)^2
W = 0.91053, p-value = 0.04192
#Not normal
##Try cubed transformation
hist((TolData_W2$Score_Full.prop)^3)
shapiro.test((TolData_W2$Score_Full.prop)^3)
Shapiro-Wilk normality test
data: (TolData_W2$Score_Full.prop)^3
W = 0.91348, p-value = 0.04837
#Not normal
##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_ScoreF_W2_lm<-lm(Score_Full.prop~Site+Genotype+Site:Genotype, data=TolData_W2)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreF_W2_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreF_W2_lm)); qqline(resid(Tol_ScoreF_W2_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_W2_lm), resid(Tol_ScoreF_W2_lm))
#Model Results
summary(Tol_ScoreF_W2_lm)
Call:
lm(formula = Score_Full.prop ~ Site + Genotype + Site:Genotype,
data = TolData_W2)
Residuals:
Min 1Q Median 3Q Max
-0.073750 -0.015096 0.002375 0.016863 0.047825
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.908519 0.007346 123.672 < 2e-16 ***
Site.L 0.019049 0.010389 1.834 0.0843 .
Genotype.L -0.074108 0.012890 -5.749 2.36e-05 ***
Genotype.Q -0.018301 0.012555 -1.458 0.1632
Site.L:Genotype.L 0.003346 0.018230 0.184 0.8565
Site.L:Genotype.Q 0.031449 0.017756 1.771 0.0945 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.03503 on 17 degrees of freedom
Multiple R-squared: 0.72, Adjusted R-squared: 0.6377
F-statistic: 8.744 on 5 and 17 DF, p-value: 0.0002968
anova(Tol_ScoreF_W2_lm)
Analysis of Variance Table
Response: Score_Full.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.005660 0.0056602 4.6129 0.04644 *
Genotype 2 0.044122 0.0220608 17.9790 6.388e-05 ***
Site:Genotype 2 0.003862 0.0019311 1.5738 0.23601
Residuals 17 0.020859 0.0012270
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreF_W2_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.08 | [0.00, 1.00]
Genotype | 0.59 | [0.28, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreF_W2_lm.res<-data.frame(anova(Tol_ScoreF_W2_lm))
Tol_ScoreF_W2_lm.res$Predictor<-rownames(Tol_ScoreF_W2_lm.res)
Tol_ScoreF_W2_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_W2_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_W2_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_W2_lm.res))
Tol_ScoreF_W2_lm.res$TimeP<-rep("W2", nrow(Tol_ScoreF_W2_lm.res))
Tol_ScoreF_W2_lm.res<-Tol_ScoreF_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreF_W2_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.933 0.0175 17 0.896 0.970
AC12 0.928 0.0175 17 0.891 0.965
AC8 0.824 0.0175 17 0.787 0.861
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.974 0.0175 17 0.937 1.011
AC12 0.919 0.0175 17 0.882 0.956
AC8 0.873 0.0202 17 0.830 0.916
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.00443 0.0248 17 0.179 0.9826
AC10 - AC8 0.10815 0.0248 17 4.366 0.0012
AC12 - AC8 0.10372 0.0248 17 4.188 0.0017
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.05555 0.0248 17 2.243 0.0923
AC10 - AC8 0.10146 0.0268 17 3.792 0.0039
AC12 - AC8 0.04591 0.0268 17 1.716 0.2281
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreF_W2_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.933 0.0175 17 0.896 0.970
SS 0.974 0.0175 17 0.937 1.011
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.928 0.0175 17 0.891 0.965
SS 0.919 0.0175 17 0.882 0.956
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.824 0.0175 17 0.787 0.861
SS 0.873 0.0202 17 0.830 0.916
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.04175 0.0248 17 -1.686 0.1101
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.00937 0.0248 17 0.378 0.7097
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.04844 0.0268 17 -1.811 0.0879
##Save p-values
#Genotypes within Sites
Tol_ScoreF_W2_lm.geno<-data.frame(emmeans(Tol_ScoreF_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreF_W2_lm.geno<-Tol_ScoreF_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_W2_lm.geno$group1<-paste(Tol_ScoreF_W2_lm.geno$Site, Tol_ScoreF_W2_lm.geno$group1, sep="_")
Tol_ScoreF_W2_lm.geno$group2<-paste(Tol_ScoreF_W2_lm.geno$Site, Tol_ScoreF_W2_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreF_W2_lm.site<-data.frame(emmeans(Tol_ScoreF_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreF_W2_lm.site<-Tol_ScoreF_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_W2_lm.site$group1<-paste(Tol_ScoreF_W2_lm.site$group1, Tol_ScoreF_W2_lm.site$Genotype, sep="_")
Tol_ScoreF_W2_lm.site$group2<-paste(Tol_ScoreF_W2_lm.site$group2, Tol_ScoreF_W2_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreF_W2_lm.p<-rbind(Tol_ScoreF_W2_lm.geno[,c(1:2,4:8)], Tol_ScoreF_W2_lm.site[,c(1:2,4:8)])
Tol_ScoreF_W2_lm.p<-Tol_ScoreF_W2_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreF_W2_lm.p$Sig<-ifelse(Tol_ScoreF_W2_lm.p$p<0.001, "***", ifelse(Tol_ScoreF_W2_lm.p$p<0.01, "**", ifelse(Tol_ScoreF_W2_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreF_W2_lm.p$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_W2_lm.p))
Tol_ScoreF_W2_lm.p$TimeP<-rep("W2", nrow(Tol_ScoreF_W2_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreF_W2_SG<-summarySE(TolData_W2, measurevar="Score_Full.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreF_W2_SG.plot<-ggplot(Tol_ScoreF_W2_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Color Retained (Full)")+
ylim(0.6, 1.2)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreF_W2_lm.p, y.position=1, step.increase=0.3, label="Sig", hide.ns=TRUE); Tol_ScoreF_W2_SG.plot
##Check normality
hist(TolData_M1$Score_Full.prop)
shapiro.test(TolData_M1$Score_Full.prop)
Shapiro-Wilk normality test
data: TolData_M1$Score_Full.prop
W = 0.9299, p-value = 0.1221
#Normal
##Model as a function of Site and Genotype
Tol_ScoreF_M1_lm<-lm(Score_Full.prop~Site+Genotype+Site:Genotype, data=TolData_M1)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreF_M1_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreF_M1_lm)); qqline(resid(Tol_ScoreF_M1_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_M1_lm), resid(Tol_ScoreF_M1_lm))
#Model Results
summary(Tol_ScoreF_M1_lm)
Call:
lm(formula = Score_Full.prop ~ Site + Genotype + Site:Genotype,
data = TolData_M1)
Residuals:
Min 1Q Median 3Q Max
-0.10493 -0.01089 0.00435 0.01777 0.11147
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.88576 0.01035 85.561 < 2e-16 ***
Site.L 0.01338 0.01464 0.914 0.374264
Genotype.L -0.07213 0.01701 -4.240 0.000623 ***
Genotype.Q 0.05226 0.01881 2.779 0.013406 *
Site.L:Genotype.L 0.02579 0.02406 1.072 0.299648
Site.L:Genotype.Q 0.02652 0.02660 0.997 0.333497
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.04811 on 16 degrees of freedom
Multiple R-squared: 0.6442, Adjusted R-squared: 0.533
F-statistic: 5.794 on 5 and 16 DF, p-value: 0.003073
anova(Tol_ScoreF_M1_lm)
Analysis of Variance Table
Response: Score_Full.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.002592 0.0025921 1.1197 0.3056920
Genotype 2 0.059505 0.0297527 12.8527 0.0004693 ***
Site:Genotype 2 0.004962 0.0024810 1.0718 0.3657507
Residuals 16 0.037038 0.0023149
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreF_M1_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.02 | [0.00, 1.00]
Genotype | 0.57 | [0.24, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreF_M1_lm.res<-data.frame(anova(Tol_ScoreF_M1_lm))
Tol_ScoreF_M1_lm.res$Predictor<-rownames(Tol_ScoreF_M1_lm.res)
Tol_ScoreF_M1_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_M1_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_M1_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M1_lm.res))
Tol_ScoreF_M1_lm.res$TimeP<-rep("M1", nrow(Tol_ScoreF_M1_lm.res))
Tol_ScoreF_M1_lm.res<-Tol_ScoreF_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreF_M1_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.954 0.0241 16 0.903 1.005
AC12 0.849 0.0278 16 0.790 0.908
AC8 0.826 0.0241 16 0.775 0.877
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.962 0.0241 16 0.911 1.013
AC12 0.837 0.0278 16 0.778 0.896
AC8 0.886 0.0241 16 0.835 0.937
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.1049 0.0367 16 2.856 0.0292
AC10 - AC8 0.1278 0.0340 16 3.756 0.0046
AC12 - AC8 0.0229 0.0367 16 0.622 0.8103
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.1251 0.0367 16 3.404 0.0096
AC10 - AC8 0.0762 0.0340 16 2.241 0.0945
AC12 - AC8 -0.0489 0.0367 16 -1.330 0.3998
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreF_M1_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.954 0.0241 16 0.903 1.005
SS 0.962 0.0241 16 0.911 1.013
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.849 0.0278 16 0.790 0.908
SS 0.837 0.0278 16 0.778 0.896
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.826 0.0241 16 0.775 0.877
SS 0.886 0.0241 16 0.835 0.937
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.00845 0.0340 16 -0.248 0.8070
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.01170 0.0393 16 0.298 0.7697
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.06003 0.0340 16 -1.764 0.0968
##Save p-values
#Genotypes within Sites
Tol_ScoreF_M1_lm.geno<-data.frame(emmeans(Tol_ScoreF_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreF_M1_lm.geno<-Tol_ScoreF_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M1_lm.geno$group1<-paste(Tol_ScoreF_M1_lm.geno$Site, Tol_ScoreF_M1_lm.geno$group1, sep="_")
Tol_ScoreF_M1_lm.geno$group2<-paste(Tol_ScoreF_M1_lm.geno$Site, Tol_ScoreF_M1_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreF_M1_lm.site<-data.frame(emmeans(Tol_ScoreF_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreF_M1_lm.site<-Tol_ScoreF_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M1_lm.site$group1<-paste(Tol_ScoreF_M1_lm.site$group1, Tol_ScoreF_M1_lm.site$Genotype, sep="_")
Tol_ScoreF_M1_lm.site$group2<-paste(Tol_ScoreF_M1_lm.site$group2, Tol_ScoreF_M1_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreF_M1_lm.p<-rbind(Tol_ScoreF_M1_lm.geno[,c(1:2,4:8)], Tol_ScoreF_M1_lm.site[,c(1:2,4:8)])
Tol_ScoreF_M1_lm.p<-Tol_ScoreF_M1_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreF_M1_lm.p$Sig<-ifelse(Tol_ScoreF_M1_lm.p$p<0.001, "***", ifelse(Tol_ScoreF_M1_lm.p$p<0.01, "**", ifelse(Tol_ScoreF_M1_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreF_M1_lm.p$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M1_lm.p))
Tol_ScoreF_M1_lm.p$TimeP<-rep("M1", nrow(Tol_ScoreF_M1_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreF_M1_SG<-summarySE(TolData_M1, measurevar="Score_Full.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreF_M1_SG.plot<-ggplot(Tol_ScoreF_M1_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Color Retained (Full)")+
ylim(0.6, 1.2)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreF_M1_lm.p, y.position=1, step.increase=0.3, label="Sig", hide.ns=TRUE); Tol_ScoreF_M1_SG.plot
##Check normality
hist(TolData_M4$Score_Full.prop)
shapiro.test(TolData_M4$Score_Full.prop)
Shapiro-Wilk normality test
data: TolData_M4$Score_Full.prop
W = 0.89096, p-value = 0.01391
#Not Normal
##Try square transformation
hist((TolData_M4$Score_Full.prop)^2)
shapiro.test((TolData_M4$Score_Full.prop)^2)
Shapiro-Wilk normality test
data: (TolData_M4$Score_Full.prop)^2
W = 0.88924, p-value = 0.01282
#Not Normal
##Try cubed transformation
hist((TolData_M4$Score_Full.prop)^3)
shapiro.test((TolData_M4$Score_Full.prop)^3)
Shapiro-Wilk normality test
data: (TolData_M4$Score_Full.prop)^3
W = 0.88739, p-value = 0.01175
#Not Normal
##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_ScoreF_M4_lm<-lm(Score_Full.prop~Site+Genotype+Site:Genotype, data=TolData_M4)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreF_M4_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreF_M4_lm)); qqline(resid(Tol_ScoreF_M4_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_M4_lm), resid(Tol_ScoreF_M4_lm))
#Model Results
summary(Tol_ScoreF_M4_lm)
Call:
lm(formula = Score_Full.prop ~ Site + Genotype + Site:Genotype,
data = TolData_M4)
Residuals:
Min 1Q Median 3Q Max
-0.058700 -0.006500 -0.000900 0.008625 0.047300
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.951088 0.004804 197.973 < 2e-16 ***
Site.L -0.001939 0.006794 -0.285 0.779
Genotype.L -0.054191 0.008321 -6.513 4.01e-06 ***
Genotype.Q 0.003016 0.008321 0.362 0.721
Site.L:Genotype.L -0.004362 0.011768 -0.371 0.715
Site.L:Genotype.Q -0.005362 0.011768 -0.456 0.654
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.02354 on 18 degrees of freedom
Multiple R-squared: 0.7048, Adjusted R-squared: 0.6228
F-statistic: 8.594 on 5 and 18 DF, p-value: 0.0002645
anova(Tol_ScoreF_M4_lm)
Analysis of Variance Table
Response: Score_Full.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.0000451 0.0000451 0.0814 0.7786
Genotype 2 0.0235660 0.0117830 21.2724 1.815e-05 ***
Site:Genotype 2 0.0001911 0.0000956 0.1725 0.8429
Residuals 18 0.0099704 0.0005539
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreF_M4_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
---------------------------------------
Site | 1.34e-03 | [0.00, 1.00]
Genotype | 0.70 | [0.45, 1.00]
Site:Genotype | 5.66e-03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreF_M4_lm.res<-data.frame(anova(Tol_ScoreF_M4_lm))
Tol_ScoreF_M4_lm.res$Predictor<-rownames(Tol_ScoreF_M4_lm.res)
Tol_ScoreF_M4_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_M4_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_M4_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M4_lm.res))
Tol_ScoreF_M4_lm.res$TimeP<-rep("M4", nrow(Tol_ScoreF_M4_lm.res))
Tol_ScoreF_M4_lm.res<-Tol_ScoreF_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreF_M4_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.991 0.0118 18 0.967 1.016
AC12 0.947 0.0118 18 0.922 0.972
AC8 0.919 0.0118 18 0.894 0.944
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.990 0.0118 18 0.965 1.015
AC12 0.950 0.0118 18 0.926 0.975
AC8 0.909 0.0118 18 0.884 0.934
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0445 0.0166 18 2.672 0.0393
AC10 - AC8 0.0723 0.0166 18 4.343 0.0011
AC12 - AC8 0.0278 0.0166 18 1.670 0.2434
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0396 0.0166 18 2.377 0.0705
AC10 - AC8 0.0810 0.0166 18 4.867 0.0003
AC12 - AC8 0.0415 0.0166 18 2.491 0.0564
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreF_M4_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.991 0.0118 18 0.967 1.016
SS 0.990 0.0118 18 0.965 1.015
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.947 0.0118 18 0.922 0.972
SS 0.950 0.0118 18 0.926 0.975
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.919 0.0118 18 0.894 0.944
SS 0.909 0.0118 18 0.884 0.934
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS 0.00147 0.0166 18 0.089 0.9304
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.00345 0.0166 18 -0.207 0.8381
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS 0.01020 0.0166 18 0.613 0.5476
##Save p-values
#Genotypes within Sites
Tol_ScoreF_M4_lm.geno<-data.frame(emmeans(Tol_ScoreF_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreF_M4_lm.geno<-Tol_ScoreF_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M4_lm.geno$group1<-paste(Tol_ScoreF_M4_lm.geno$Site, Tol_ScoreF_M4_lm.geno$group1, sep="_")
Tol_ScoreF_M4_lm.geno$group2<-paste(Tol_ScoreF_M4_lm.geno$Site, Tol_ScoreF_M4_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreF_M4_lm.site<-data.frame(emmeans(Tol_ScoreF_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreF_M4_lm.site<-Tol_ScoreF_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M4_lm.site$group1<-paste(Tol_ScoreF_M4_lm.site$group1, Tol_ScoreF_M4_lm.site$Genotype, sep="_")
Tol_ScoreF_M4_lm.site$group2<-paste(Tol_ScoreF_M4_lm.site$group2, Tol_ScoreF_M4_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreF_M4_lm.p<-rbind(Tol_ScoreF_M4_lm.geno[,c(1:2,4:8)], Tol_ScoreF_M4_lm.site[,c(1:2,4:8)])
Tol_ScoreF_M4_lm.p<-Tol_ScoreF_M4_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreF_M4_lm.p$Sig<-ifelse(Tol_ScoreF_M4_lm.p$p<0.001, "***", ifelse(Tol_ScoreF_M4_lm.p$p<0.01, "**", ifelse(Tol_ScoreF_M4_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreF_M4_lm.p$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M4_lm.p))
Tol_ScoreF_M4_lm.p$TimeP<-rep("M4", nrow(Tol_ScoreF_M4_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreF_M4_SG<-summarySE(TolData_M4, measurevar="Score_Full.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreF_M4_SG.plot<-ggplot(Tol_ScoreF_M4_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Color Retained (Full)")+
ylim(0.6, 1.2)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreF_M4_lm.p, y.position=1.05, step.increase=0.35, label="Sig", hide.ns=TRUE); Tol_ScoreF_M4_SG.plot
##Check normality
hist(TolData$Score_TP.prop)
shapiro.test(TolData$Score_TP.prop)
Shapiro-Wilk normality test
data: TolData$Score_TP.prop
W = 0.95179, p-value = 0.009694
#Not Normal
##Try square transformation
hist((TolData$Score_TP.prop)^2)
shapiro.test((TolData$Score_TP.prop)^2)
Shapiro-Wilk normality test
data: (TolData$Score_TP.prop)^2
W = 0.95179, p-value = 0.009688
#Not normal
##Try cubed transformation
hist((TolData$Score_TP.prop)^3)
shapiro.test((TolData$Score_TP.prop)^3)
Shapiro-Wilk normality test
data: (TolData$Score_TP.prop)^3
W = 0.95056, p-value = 0.00836
#Not normal
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreTP_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreTP_lm)); qqline(resid(Tol_ScoreTP_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_lm), resid(Tol_ScoreTP_lm))
#Model Results
summary(Tol_ScoreTP_lm)
Call:
lm(formula = Score_TP.prop ~ Site + Genotype + TimeP + Site:Genotype +
Site:TimeP + Genotype:TimeP, data = TolData)
Residuals:
Min 1Q Median 3Q Max
-0.133294 -0.019060 0.000725 0.026225 0.122906
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.883494 0.006319 139.809 < 2e-16 ***
Site.L 0.038546 0.008912 4.325 6.47e-05 ***
Genotype.L -0.102221 0.010797 -9.468 3.81e-13 ***
Genotype.Q 0.004971 0.011092 0.448 0.655784
TimeP.L 0.068495 0.010797 6.344 4.45e-08 ***
TimeP.Q 0.021265 0.011092 1.917 0.060411 .
Site.L:Genotype.L 0.014004 0.015269 0.917 0.363050
Site.L:Genotype.Q -0.008574 0.015625 -0.549 0.585424
Site.L:TimeP.L -0.036792 0.015269 -2.410 0.019344 *
Site.L:TimeP.Q -0.004886 0.015625 -0.313 0.755697
Genotype.L:TimeP.L 0.029281 0.018818 1.556 0.125448
Genotype.Q:TimeP.L 0.024642 0.018582 1.326 0.190292
Genotype.L:TimeP.Q 0.015340 0.018582 0.826 0.412653
Genotype.Q:TimeP.Q -0.072757 0.019821 -3.671 0.000548 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05222 on 55 degrees of freedom
Multiple R-squared: 0.763, Adjusted R-squared: 0.7069
F-statistic: 13.62 on 13 and 55 DF, p-value: 8.098e-13
anova(Tol_ScoreTP_lm)
Analysis of Variance Table
Response: Score_TP.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.057199 0.057199 20.9741 2.709e-05 ***
Genotype 2 0.238412 0.119206 43.7110 4.331e-12 ***
TimeP 2 0.117103 0.058552 21.4700 1.284e-07 ***
Site:Genotype 2 0.003221 0.001611 0.5906 0.557481
Site:TimeP 2 0.017283 0.008642 3.1687 0.049832 *
Genotype:TimeP 4 0.049580 0.012395 4.5450 0.003041 **
Residuals 55 0.149993 0.002727
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreTP_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
----------------------------------------
Site | 0.09 | [0.01, 1.00]
Genotype | 0.38 | [0.20, 1.00]
TimeP | 0.19 | [0.04, 1.00]
Site:Genotype | 5.09e-03 | [0.00, 1.00]
Site:TimeP | 0.03 | [0.00, 1.00]
Genotype:TimeP | 0.08 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreTP_lm.res<-data.frame(anova(Tol_ScoreTP_lm))
Tol_ScoreTP_lm.res$Predictor<-rownames(Tol_ScoreTP_lm.res)
Tol_ScoreTP_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_lm.res))
Tol_ScoreTP_lm.res<-Tol_ScoreTP_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Strong influence of Timepoint, including interactions with variables of interest, Site and Genotype. Will analyze each Timepoint individually.
##Check normality
hist(TolData_W2$Score_TP.prop)
shapiro.test(TolData_W2$Score_TP.prop)
Shapiro-Wilk normality test
data: TolData_W2$Score_TP.prop
W = 0.95707, p-value = 0.4068
#Normal
##Model as a function of Site and Genotype
Tol_ScoreTP_W2_lm<-lm(Score_TP.prop~Site+Genotype+Site:Genotype, data=TolData_W2)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreTP_W2_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreTP_W2_lm)); qqline(resid(Tol_ScoreTP_W2_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_W2_lm), resid(Tol_ScoreTP_W2_lm))
#Model Results
summary(Tol_ScoreTP_W2_lm)
Call:
lm(formula = Score_TP.prop ~ Site + Genotype + Site:Genotype,
data = TolData_W2)
Residuals:
Min 1Q Median 3Q Max
-0.06045 -0.02293 0.00115 0.02314 0.07585
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.8997236 0.0084985 105.868 < 2e-16 ***
Site.L 0.0162968 0.0120187 1.356 0.192853
Genotype.L -0.0664798 0.0149123 -4.458 0.000345 ***
Genotype.Q -0.0201334 0.0145249 -1.386 0.183620
Site.L:Genotype.L 0.0005833 0.0210892 0.028 0.978255
Site.L:Genotype.Q 0.0026173 0.0205413 0.127 0.900105
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.04052 on 17 degrees of freedom
Multiple R-squared: 0.5884, Adjusted R-squared: 0.4674
F-statistic: 4.861 on 5 and 17 DF, p-value: 0.006075
anova(Tol_ScoreTP_W2_lm)
Analysis of Variance Table
Response: Score_TP.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.004491 0.0044907 2.7346 0.1165399
Genotype 2 0.035398 0.0176991 10.7779 0.0009485 ***
Site:Genotype 2 0.000027 0.0000137 0.0084 0.9916821
Residuals 17 0.027917 0.0016422
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreTP_W2_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
---------------------------------------
Site | 0.07 | [0.00, 1.00]
Genotype | 0.52 | [0.19, 1.00]
Site:Genotype | 4.05e-04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreTP_W2_lm.res<-data.frame(anova(Tol_ScoreTP_W2_lm))
Tol_ScoreTP_W2_lm.res$Predictor<-rownames(Tol_ScoreTP_W2_lm.res)
Tol_ScoreTP_W2_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_W2_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_W2_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_W2_lm.res))
Tol_ScoreTP_W2_lm.res$TimeP<-rep("W2", nrow(Tol_ScoreTP_W2_lm.res))
Tol_ScoreTP_W2_lm.res<-Tol_ScoreTP_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreTP_W2_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.927 0.0203 17 0.884 0.969
AC12 0.906 0.0203 17 0.863 0.949
AC8 0.832 0.0203 17 0.789 0.875
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.951 0.0203 17 0.908 0.993
AC12 0.926 0.0203 17 0.883 0.969
AC8 0.857 0.0234 17 0.808 0.906
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0204 0.0287 17 0.711 0.7604
AC10 - AC8 0.0946 0.0287 17 3.301 0.0111
AC12 - AC8 0.0742 0.0287 17 2.590 0.0476
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0243 0.0287 17 0.849 0.6786
AC10 - AC8 0.0934 0.0310 17 3.019 0.0201
AC12 - AC8 0.0691 0.0310 17 2.233 0.0940
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreTP_W2_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.927 0.0203 17 0.884 0.969
SS 0.951 0.0203 17 0.908 0.993
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.906 0.0203 17 0.863 0.949
SS 0.926 0.0203 17 0.883 0.969
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.832 0.0203 17 0.789 0.875
SS 0.857 0.0234 17 0.808 0.906
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0240 0.0287 17 -0.837 0.4144
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.0200 0.0287 17 -0.699 0.4941
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.0251 0.0310 17 -0.812 0.4278
##Save p-values
#Genotypes within Sites
Tol_ScoreTP_W2_lm.geno<-data.frame(emmeans(Tol_ScoreTP_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreTP_W2_lm.geno<-Tol_ScoreTP_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_W2_lm.geno$group1<-paste(Tol_ScoreTP_W2_lm.geno$Site, Tol_ScoreTP_W2_lm.geno$group1, sep="_")
Tol_ScoreTP_W2_lm.geno$group2<-paste(Tol_ScoreTP_W2_lm.geno$Site, Tol_ScoreTP_W2_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreTP_W2_lm.site<-data.frame(emmeans(Tol_ScoreTP_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreTP_W2_lm.site<-Tol_ScoreTP_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_W2_lm.site$group1<-paste(Tol_ScoreTP_W2_lm.site$group1, Tol_ScoreTP_W2_lm.site$Genotype, sep="_")
Tol_ScoreTP_W2_lm.site$group2<-paste(Tol_ScoreTP_W2_lm.site$group2, Tol_ScoreTP_W2_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreTP_W2_lm.p<-rbind(Tol_ScoreTP_W2_lm.geno[,c(1:2,4:8)], Tol_ScoreTP_W2_lm.site[,c(1:2,4:8)])
Tol_ScoreTP_W2_lm.p<-Tol_ScoreTP_W2_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreTP_W2_lm.p$Sig<-ifelse(Tol_ScoreTP_W2_lm.p$p<0.001, "***", ifelse(Tol_ScoreTP_W2_lm.p$p<0.01, "**", ifelse(Tol_ScoreTP_W2_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreTP_W2_lm.p$Response<-rep("Color_TP", nrow(Tol_ScoreTP_W2_lm.p))
Tol_ScoreTP_W2_lm.p$TimeP<-rep("W2", nrow(Tol_ScoreTP_W2_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreTP_W2_SG<-summarySE(TolData_W2, measurevar="Score_TP.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreTP_W2_SG.plot<-ggplot(Tol_ScoreTP_W2_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Color Retained (TP)")+
ylim(0.6, 1.2)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreTP_W2_lm.p, y.position=1, step.increase=0.3, label="Sig", hide.ns=TRUE); Tol_ScoreTP_W2_SG.plot
##Check normality
hist(TolData_M1$Score_TP.prop)
shapiro.test(TolData_M1$Score_TP.prop)
Shapiro-Wilk normality test
data: TolData_M1$Score_TP.prop
W = 0.92691, p-value = 0.1058
#Normal
##Model as a function of Site and Genotype
Tol_ScoreTP_M1_lm<-lm(Score_TP.prop~Site+Genotype+Site:Genotype, data=TolData_M1)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreTP_M1_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreTP_M1_lm)); qqline(resid(Tol_ScoreTP_M1_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_M1_lm), resid(Tol_ScoreTP_M1_lm))
#Model Results
summary(Tol_ScoreTP_M1_lm)
Call:
lm(formula = Score_TP.prop ~ Site + Genotype + Site:Genotype,
data = TolData_M1)
Residuals:
Min 1Q Median 3Q Max
-0.049233 -0.015825 0.003625 0.011990 0.071567
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.880342 0.006598 133.425 < 2e-16 ***
Site.L 0.020003 0.009331 2.144 0.04776 *
Genotype.L -0.076173 0.010842 -7.026 2.86e-06 ***
Genotype.Q 0.047387 0.011986 3.954 0.00114 **
Site.L:Genotype.L 0.019850 0.015332 1.295 0.21382
Site.L:Genotype.Q 0.025394 0.016951 1.498 0.15357
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.03066 on 16 degrees of freedom
Multiple R-squared: 0.8232, Adjusted R-squared: 0.768
F-statistic: 14.9 on 5 and 16 DF, p-value: 1.559e-05
anova(Tol_ScoreTP_M1_lm)
Analysis of Variance Table
Response: Score_TP.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.005270 0.0052700 5.6044 0.03086 *
Genotype 2 0.061117 0.0305585 32.4976 2.319e-06 ***
Site:Genotype 2 0.003686 0.0018432 1.9602 0.17321
Residuals 16 0.015045 0.0009403
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreTP_M1_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.06 | [0.00, 1.00]
Genotype | 0.72 | [0.47, 1.00]
Site:Genotype | 0.04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreTP_M1_lm.res<-data.frame(anova(Tol_ScoreTP_M1_lm))
Tol_ScoreTP_M1_lm.res$Predictor<-rownames(Tol_ScoreTP_M1_lm.res)
Tol_ScoreTP_M1_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_M1_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_M1_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M1_lm.res))
Tol_ScoreTP_M1_lm.res$TimeP<-rep("M1", nrow(Tol_ScoreTP_M1_lm.res))
Tol_ScoreTP_M1_lm.res<-Tol_ScoreTP_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreTP_M1_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.942 0.0153 16 0.909 0.975
AC12 0.842 0.0177 16 0.805 0.880
AC8 0.814 0.0153 16 0.782 0.847
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.965 0.0153 16 0.933 0.998
AC12 0.841 0.0177 16 0.804 0.879
AC8 0.877 0.0153 16 0.845 0.910
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0998 0.0234 16 4.263 0.0016
AC10 - AC8 0.1276 0.0217 16 5.884 0.0001
AC12 - AC8 0.0277 0.0234 16 1.184 0.4789
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.1240 0.0234 16 5.293 0.0002
AC10 - AC8 0.0879 0.0217 16 4.053 0.0025
AC12 - AC8 -0.0361 0.0234 16 -1.541 0.2990
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreTP_M1_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.942 0.0153 16 0.909 0.975
SS 0.965 0.0153 16 0.933 0.998
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.842 0.0177 16 0.805 0.880
SS 0.841 0.0177 16 0.804 0.879
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.814 0.0153 16 0.782 0.847
SS 0.877 0.0153 16 0.845 0.910
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.02310 0.0217 16 -1.065 0.3025
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.00103 0.0250 16 0.041 0.9676
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.06280 0.0217 16 -2.896 0.0105
##Save p-values
#Genotypes within Sites
Tol_ScoreTP_M1_lm.geno<-data.frame(emmeans(Tol_ScoreTP_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreTP_M1_lm.geno<-Tol_ScoreTP_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M1_lm.geno$group1<-paste(Tol_ScoreTP_M1_lm.geno$Site, Tol_ScoreTP_M1_lm.geno$group1, sep="_")
Tol_ScoreTP_M1_lm.geno$group2<-paste(Tol_ScoreTP_M1_lm.geno$Site, Tol_ScoreTP_M1_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreTP_M1_lm.site<-data.frame(emmeans(Tol_ScoreTP_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreTP_M1_lm.site<-Tol_ScoreTP_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M1_lm.site$group1<-paste(Tol_ScoreTP_M1_lm.site$group1, Tol_ScoreTP_M1_lm.site$Genotype, sep="_")
Tol_ScoreTP_M1_lm.site$group2<-paste(Tol_ScoreTP_M1_lm.site$group2, Tol_ScoreTP_M1_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreTP_M1_lm.p<-rbind(Tol_ScoreTP_M1_lm.geno[,c(1:2,4:8)], Tol_ScoreTP_M1_lm.site[,c(1:2,4:8)])
Tol_ScoreTP_M1_lm.p<-Tol_ScoreTP_M1_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreTP_M1_lm.p$Sig<-ifelse(Tol_ScoreTP_M1_lm.p$p<0.001, "***", ifelse(Tol_ScoreTP_M1_lm.p$p<0.01, "**", ifelse(Tol_ScoreTP_M1_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreTP_M1_lm.p$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M1_lm.p))
Tol_ScoreTP_M1_lm.p$TimeP<-rep("M1", nrow(Tol_ScoreTP_M1_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreTP_M1_SG<-summarySE(TolData_M1, measurevar="Score_TP.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreTP_M1_SG.plot<-ggplot(Tol_ScoreTP_M1_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Color Retained (TP)")+
ylim(0.6, 1.2)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreTP_M1_lm.p, y.position=1, step.increase=0.25, label="Sig", hide.ns=TRUE); Tol_ScoreTP_M1_SG.plot
##Check normality
hist(TolData_M4$Score_TP.prop)
shapiro.test(TolData_M4$Score_TP.prop)
Shapiro-Wilk normality test
data: TolData_M4$Score_TP.prop
W = 0.92944, p-value = 0.09471
#Normal
##Model as a function of Site and Genotype
Tol_ScoreTP_M4_lm<-lm(Score_TP.prop~Site+Genotype+Site:Genotype, data=TolData_M4)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreTP_M4_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreTP_M4_lm)); qqline(resid(Tol_ScoreTP_M4_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_M4_lm), resid(Tol_ScoreTP_M4_lm))
#Model Results
summary(Tol_ScoreTP_M4_lm)
Call:
lm(formula = Score_TP.prop ~ Site + Genotype + Site:Genotype,
data = TolData_M4)
Residuals:
Min 1Q Median 3Q Max
-0.107100 -0.019856 0.004288 0.021244 0.073300
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.914738 0.009447 96.827 < 2e-16 ***
Site.L 0.002646 0.013360 0.198 0.845240
Genotype.L -0.073451 0.016363 -4.489 0.000284 ***
Genotype.Q 0.001990 0.016363 0.122 0.904540
Site.L:Genotype.L -0.016375 0.023141 -0.708 0.488238
Site.L:Genotype.Q 0.026775 0.023141 1.157 0.262382
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.04628 on 18 degrees of freedom
Multiple R-squared: 0.5505, Adjusted R-squared: 0.4256
F-statistic: 4.409 on 5 and 18 DF, p-value: 0.008508
anova(Tol_ScoreTP_M4_lm)
Analysis of Variance Table
Response: Score_TP.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.000084 0.000084 0.0392 0.845240
Genotype 2 0.043192 0.021596 10.0823 0.001155 **
Site:Genotype 2 0.003940 0.001970 0.9197 0.416564
Residuals 18 0.038555 0.002142
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreTP_M4_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
---------------------------------------
Site | 9.79e-04 | [0.00, 1.00]
Genotype | 0.50 | [0.18, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreTP_M4_lm.res<-data.frame(anova(Tol_ScoreTP_M4_lm))
Tol_ScoreTP_M4_lm.res$Predictor<-rownames(Tol_ScoreTP_M4_lm.res)
Tol_ScoreTP_M4_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_M4_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_M4_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M4_lm.res))
Tol_ScoreTP_M4_lm.res$TimeP<-rep("M4", nrow(Tol_ScoreTP_M4_lm.res))
Tol_ScoreTP_M4_lm.res<-Tol_ScoreTP_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreTP_M4_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.950 0.0231 18 0.901 0.998
AC12 0.927 0.0231 18 0.878 0.975
AC8 0.862 0.0231 18 0.814 0.911
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.985 0.0231 18 0.937 1.034
AC12 0.900 0.0231 18 0.851 0.948
AC8 0.865 0.0231 18 0.816 0.914
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0230 0.0327 18 0.703 0.7649
AC10 - AC8 0.0875 0.0327 18 2.674 0.0392
AC12 - AC8 0.0645 0.0327 18 1.971 0.1482
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0858 0.0327 18 2.620 0.0437
AC10 - AC8 0.1202 0.0327 18 3.674 0.0047
AC12 - AC8 0.0345 0.0327 18 1.054 0.5536
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreTP_M4_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.950 0.0231 18 0.901 0.998
SS 0.985 0.0231 18 0.937 1.034
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.927 0.0231 18 0.878 0.975
SS 0.900 0.0231 18 0.851 0.948
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.862 0.0231 18 0.814 0.911
SS 0.865 0.0231 18 0.816 0.914
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.03558 0.0327 18 -1.087 0.2914
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.02718 0.0327 18 0.830 0.4172
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.00282 0.0327 18 -0.086 0.9322
##Save p-values
#Genotypes within Sites
Tol_ScoreTP_M4_lm.geno<-data.frame(emmeans(Tol_ScoreTP_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreTP_M4_lm.geno<-Tol_ScoreTP_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M4_lm.geno$group1<-paste(Tol_ScoreTP_M4_lm.geno$Site, Tol_ScoreTP_M4_lm.geno$group1, sep="_")
Tol_ScoreTP_M4_lm.geno$group2<-paste(Tol_ScoreTP_M4_lm.geno$Site, Tol_ScoreTP_M4_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreTP_M4_lm.site<-data.frame(emmeans(Tol_ScoreTP_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreTP_M4_lm.site<-Tol_ScoreTP_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M4_lm.site$group1<-paste(Tol_ScoreTP_M4_lm.site$group1, Tol_ScoreTP_M4_lm.site$Genotype, sep="_")
Tol_ScoreTP_M4_lm.site$group2<-paste(Tol_ScoreTP_M4_lm.site$group2, Tol_ScoreTP_M4_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreTP_M4_lm.p<-rbind(Tol_ScoreTP_M4_lm.geno[,c(1:2,4:8)], Tol_ScoreTP_M4_lm.site[,c(1:2,4:8)])
Tol_ScoreTP_M4_lm.p<-Tol_ScoreTP_M4_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreTP_M4_lm.p$Sig<-ifelse(Tol_ScoreTP_M4_lm.p$p<0.001, "***", ifelse(Tol_ScoreTP_M4_lm.p$p<0.01, "**", ifelse(Tol_ScoreTP_M4_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreTP_M4_lm.p$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M4_lm.p))
Tol_ScoreTP_M4_lm.p$TimeP<-rep("M4", nrow(Tol_ScoreTP_M4_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreTP_M4_SG<-summarySE(TolData_M4, measurevar="Score_TP.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreTP_M4_SG.plot<-ggplot(Tol_ScoreTP_M4_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Color Retained (TP)")+
ylim(0.6, 1.2)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreTP_M4_lm.p, y.position=1, step.increase=0.3, label="Sig", hide.ns=TRUE); Tol_ScoreTP_M4_SG.plot
Creating a heatmap to compare the direction and significance of pairwise comparisons across Tolerance metrics. Positive estimates indicate that Group 1 > Group 2. P values of less than 0.05 are considered significant.
##Combine Results
Tol_contrasts<-rbind(Tol_Chl_W2_lm.p, Tol_Chl_M1_lm.p, Tol_Chl_M4_lm.p,
Tol_Sym_W2_lm.p, Tol_Sym_M1_lm.p, Tol_Sym_M4_lm.p,
Tol_ScoreF_W2_lm.p, Tol_ScoreF_M1_lm.p, Tol_ScoreF_M4_lm.p,
Tol_ScoreTP_W2_lm.p, Tol_ScoreTP_M1_lm.p, Tol_ScoreTP_M4_lm.p)
##Create Contrast Column
Tol_contrasts$Contrast<-paste(Tol_contrasts$group1, Tol_contrasts$group2, sep=" v. ")
##Add Set column with Contrast and Timepoint
Tol_contrasts$Set<-paste(Tol_contrasts$TimeP, Tol_contrasts$Contrast, sep=" ")
##Re-order by Seasonal Timepoint
Tol_contrasts$TimeP<-factor(Tol_contrasts$TimeP, levels=c("W2", "M1", "M4"), ordered=TRUE)
Tol_contrasts<- Tol_contrasts %>% arrange(desc(as.numeric(TimeP)))
Tol_contrasts$Set<-factor(Tol_contrasts$Set, levels=c(unique(Tol_contrasts$Set)), ordered=TRUE)
Convert Estimate to the Response scale (instead of log +1 scale) for models where the response was log transformed
Tol_contrasts$Estimate<-Tol_contrasts$estimate
#Chl M4
Tol_contrasts$Estimate[which(Tol_contrasts$Response=="Chlorophyll" & Tol_contrasts$TimeP=="M4")]<-c(exp(Tol_contrasts$estimate[which(Tol_contrasts$Response=="Chlorophyll" & Tol_contrasts$TimeP=="M4")])-1)
#Sym M1
Tol_contrasts$Estimate[which(Tol_contrasts$Response=="Symbionts" & Tol_contrasts$TimeP=="M1")]<-c(exp(Tol_contrasts$estimate[which(Tol_contrasts$Response=="Symbionts" & Tol_contrasts$TimeP=="M1")])-1)
Tol_Pairs_W2.plot<-ggplot(data=Tol_contrasts[which(Tol_contrasts$TimeP=="W2"),], aes(x=Response, y=Contrast, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Pairwise Contrasts Summer 1")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs_W2.plot
Tol_Pairs_M1.plot<-ggplot(data=Tol_contrasts[which(Tol_contrasts$TimeP=="M1"),], aes(x=Response, y=Contrast, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Pairwise Contrasts Summer 2")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs_M1.plot
Tol_Pairs_M4.plot<-ggplot(data=Tol_contrasts[which(Tol_contrasts$TimeP=="M4"),], aes(x=Response, y=Contrast, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Pairwise Contrasts Month 4")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs_M4.plot
Tol_Pairs.plot<-ggplot(data=Tol_contrasts, aes(x=Response, y=Set, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Pairwise Contrasts")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs.plot
Subsetting Contrasts to only include contrasts with significant differences in at least one of the thermal tolerance metrics (Chlorophyll or Symbiont or Color retention)
##Remove rows with non-significant p-values
Tol_contrasts_sig<-subset(Tol_contrasts, p < 0.05)
##Keep Sets where at least one metric shows a significant difference
Tol_Sets_sig<-data.frame(Set=c(unique(Tol_contrasts_sig$Set)))
Tol_contrasts_sig_Set<-merge(Tol_Sets_sig, Tol_contrasts, all.x=TRUE, all.y=FALSE)
Tol_Pairs_sig.plot<-ggplot(data=Tol_contrasts_sig_Set, aes(x=Response, y=Set, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Significant Pairwise Contrasts")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs_sig.plot
Assessing the proportion of agreement between methods as the pairwise comparisons where two methods both find a significant difference out of the total pairwise comparisons. Calculating the proportion of agreement between Chlorophyll and Color Full Set, Chlorophyll and Color Timepoint, Symbionts and Color Full Set, Symbionts and Color Timepoint, and Chlorophyll and Symbionts. Comparing the proportion of agreement between methods of quantifying thermal tolerance.
names(Tol_contrasts)
[1] "group1" "group2" "estimate" "SE" "df" "t.ratio" "p" "Sig"
[9] "Response" "TimeP" "Contrast" "Set" "Estimate"
##Change long to short format
Tol_contrasts.agree<-Tol_contrasts %>% pivot_wider(names_from="Response", values_from=c("p", "Estimate"), id_cols=c("TimeP", "Contrast", "Set"))
##Convert p values to binary
#0: No significant difference (p>0.05) detected in that pairwise contrast
#1: Significant difference (p<0.05) detected in that pairwise contrast
Tol_contrasts.agree[,c(4:7)]<-ifelse(Tol_contrasts.agree[,c(4:7)]<0.05, 1, 0)
##Convert estimates to Greater or Less
#Greater: Positive estimates indicate A>B in the Contrast
#Less: Negative estimates indicate A<B in the Contrast
Tol_contrasts.agree[,c(8:11)]<-ifelse(Tol_contrasts.agree[,c(8:11)]<0, "Less", "Greater")
##Identify agreement
#1: Both metrics identified a significant difference in that pairwise contrast OR Neither metric identified a significant difference
#0: One metrics identified a significant difference in that pairwise contrast
#For each matching case (1), confirm matching Estimates
Tol_contrasts.agree$Chl.Sym<-ifelse(Tol_contrasts.agree$p_Chlorophyll==1 &
Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Symbionts &
Tol_contrasts.agree$Estimate_Chlorophyll==Tol_contrasts.agree$Estimate_Symbionts, 1,
ifelse(Tol_contrasts.agree$p_Chlorophyll==0 &
Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Symbionts, 1, 0))
Tol_contrasts.agree$Chl.ColF<-ifelse(Tol_contrasts.agree$p_Chlorophyll==1 & Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Color_FullSet &
Tol_contrasts.agree$Estimate_Chlorophyll==Tol_contrasts.agree$Estimate_Color_FullSet, 1,
ifelse(Tol_contrasts.agree$p_Chlorophyll==0 &
Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Color_FullSet, 1, 0))
Tol_contrasts.agree$Chl.ColTP<-ifelse(Tol_contrasts.agree$p_Chlorophyll==1 &
Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Color_TP & Tol_contrasts.agree$Estimate_Chlorophyll==Tol_contrasts.agree$Estimate_Color_TP, 1,
ifelse(Tol_contrasts.agree$p_Chlorophyll==0 &
Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Color_TP, 1, 0))
Tol_contrasts.agree$Sym.ColF<-ifelse(Tol_contrasts.agree$p_Symbionts==1 &
Tol_contrasts.agree$p_Symbionts==Tol_contrasts.agree$p_Color_FullSet & Tol_contrasts.agree$Estimate_Symbionts==Tol_contrasts.agree$Estimate_Color_FullSet, 1,
ifelse(Tol_contrasts.agree$p_Symbionts==0 &
Tol_contrasts.agree$p_Symbionts==Tol_contrasts.agree$p_Color_FullSet, 1, 0))
Tol_contrasts.agree$Sym.ColTP<-ifelse(Tol_contrasts.agree$p_Symbionts==1 &
Tol_contrasts.agree$p_Symbionts==Tol_contrasts.agree$p_Color_TP & Tol_contrasts.agree$Estimate_Symbionts==Tol_contrasts.agree$Estimate_Color_TP, 1,
ifelse(Tol_contrasts.agree$p_Symbionts==0 &
Tol_contrasts.agree$p_Symbionts==Tol_contrasts.agree$p_Color_TP, 1, 0))
Tol_contrasts.agree$ColF.ColTP<-ifelse(Tol_contrasts.agree$p_Color_FullSet==1 &
Tol_contrasts.agree$p_Color_FullSet==Tol_contrasts.agree$p_Color_TP & Tol_contrasts.agree$Estimate_Color_FullSet==Tol_contrasts.agree$Estimate_Color_TP, 1,
ifelse(Tol_contrasts.agree$p_Color_FullSet==0 &
Tol_contrasts.agree$p_Color_FullSet==Tol_contrasts.agree$p_Color_TP, 1, 0))
##Convert to long to calculate instances of agreement
Tol_contrasts.agree.long<-Tol_contrasts.agree %>% pivot_longer(cols=c("Chl.Sym", "Chl.ColF", "Chl.ColTP", "Sym.ColF", "Sym.ColTP", "ColF.ColTP"), names_to="Compare", values_to="Agree")
##Sum Cases of Agreement
Tol_contrasts.agree.sum<-aggregate(Tol_contrasts.agree.long$Agree, list(Tol_contrasts.agree.long$Compare), sum)
names(Tol_contrasts.agree.sum)<-c("Compare", "Agree")
##Total Pairwise Comparisons
Tol_contrasts.agree.sum$Total<-aggregate(Tol_contrasts.agree.long$Agree, list(Tol_contrasts.agree.long$Compare), function(x) length(x))[,2]
prop.test(Tol_contrasts.agree.sum$Agree, Tol_contrasts.agree.sum$Total)
6-sample test for equality of proportions without continuity correction
data: Tol_contrasts.agree.sum$Agree out of Tol_contrasts.agree.sum$Total
X-squared = 2.8177, df = 5, p-value = 0.7281
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3 prop 4 prop 5 prop 6
0.7037037 0.7407407 0.7037037 0.8518519 0.6666667 0.7407407
No significant difference in the proportion of agreement in detecting pairwise differences in thermal tolerance across comparisons of traditional bleaching or color score methods (Pearson’s chi-squared test p-value = 0.7281)
##Sum Cases of Agreement by Comparison and Timepoint
Tol_contrasts.agree.sum.tp<-aggregate(Tol_contrasts.agree.long$Agree, list(Tol_contrasts.agree.long$TimeP, Tol_contrasts.agree.long$Compare), sum)
names(Tol_contrasts.agree.sum.tp)<-c("TimeP", "Compare", "Agree")
##Total Pairwise Comparisons
Tol_contrasts.agree.sum.tp$Total<-aggregate(Tol_contrasts.agree.long$Agree, list(Tol_contrasts.agree.long$TimeP, Tol_contrasts.agree.long$Compare), function(x) length(x))[,3]
##Subset by Comparison
Tol_contrasts.agree.Chl.Sym<-subset(Tol_contrasts.agree.sum.tp, Compare=="Chl.Sym")
##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Chl.Sym$Agree, Tol_contrasts.agree.Chl.Sym$Total)
Warning: Chi-squared approximation may be incorrect
3-sample test for equality of proportions without continuity correction
data: Tol_contrasts.agree.Chl.Sym$Agree out of Tol_contrasts.agree.Chl.Sym$Total
X-squared = 5.6842, df = 2, p-value = 0.0583
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3
0.5555556 1.0000000 0.5555556
##Subset by Comparison
Tol_contrasts.agree.Chl.ColF<-subset(Tol_contrasts.agree.sum.tp, Compare=="Chl.ColF")
##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Chl.ColF$Agree, Tol_contrasts.agree.Chl.ColF$Total)
Warning: Chi-squared approximation may be incorrect
3-sample test for equality of proportions without continuity correction
data: Tol_contrasts.agree.Chl.ColF$Agree out of Tol_contrasts.agree.Chl.ColF$Total
X-squared = 9.9474, df = 2, p-value = 0.006918
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3
1.0000000 0.7777778 0.3333333
##Subset by Comparison
Tol_contrasts.agree.Chl.ColTP<-subset(Tol_contrasts.agree.sum.tp, Compare=="Chl.ColTP")
##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Chl.ColTP$Agree, Tol_contrasts.agree.Chl.ColTP$Total)
Warning: Chi-squared approximation may be incorrect
3-sample test for equality of proportions without continuity correction
data: Tol_contrasts.agree.Chl.ColTP$Agree out of Tol_contrasts.agree.Chl.ColTP$Total
X-squared = 18.9, df = 2, p-value = 7.869e-05
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3
1.0000000 1.0000000 0.2222222
##Subset by Comparison
Tol_contrasts.agree.Sym.ColF<-subset(Tol_contrasts.agree.sum.tp, Compare=="Sym.ColF")
##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Sym.ColF$Agree, Tol_contrasts.agree.Sym.ColF$Total)
Warning: Chi-squared approximation may be incorrect
3-sample test for equality of proportions without continuity correction
data: Tol_contrasts.agree.Sym.ColF$Agree out of Tol_contrasts.agree.Sym.ColF$Total
X-squared = 1, df = 2, p-value = 0.6065
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3
0.5555556 0.7777778 0.6666667
##Subset by Comparison
Tol_contrasts.agree.Sym.ColTP<-subset(Tol_contrasts.agree.sum.tp, Compare=="Sym.ColTP")
##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Sym.ColTP$Agree, Tol_contrasts.agree.Sym.ColTP$Total)
Warning: Chi-squared approximation may be incorrect
3-sample test for equality of proportions without continuity correction
data: Tol_contrasts.agree.Sym.ColTP$Agree out of Tol_contrasts.agree.Sym.ColTP$Total
X-squared = 5.0143, df = 2, p-value = 0.0815
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3
0.5555556 1.0000000 0.6666667
##Subset by Comparison
Tol_contrasts.agree.ColF.ColTP<-subset(Tol_contrasts.agree.sum.tp, Compare=="ColF.ColTP")
##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.ColF.ColTP$Agree, Tol_contrasts.agree.ColF.ColTP$Total)
Warning: Chi-squared approximation may be incorrect
3-sample test for equality of proportions without continuity correction
data: Tol_contrasts.agree.ColF.ColTP$Agree out of Tol_contrasts.agree.ColF.ColTP$Total
X-squared = 2.3478, df = 2, p-value = 0.3092
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3
1.0000000 0.7777778 0.7777778
##Top row
Tol_Chl_W2_SG.plot<-Tol_Chl_W2_SG.plot + ggtitle("Summer 1") + labs(x="")+
theme(legend.position="top", legend.direction="horizontal", plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))
Tol_Chl_M1_SG.plot<-Tol_Chl_M1_SG.plot + ggtitle("Summer 2") + labs(x="", y="")+
theme(legend.position="top", legend.direction="horizontal", plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))
Tol_Chl_M4_SG.plot<-Tol_Chl_M4_SG.plot + ggtitle("Winter") + labs(x="", y="")+
theme(legend.position="top", legend.direction="horizontal", plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))
##Bottom row
Tol_Sym_W2_SG.plot<-Tol_Sym_W2_SG.plot+labs(x="Site and Genotype")+ theme(legend.position="none")
Tol_Sym_M1_SG.plot<-Tol_Sym_M1_SG.plot+labs(x="Site and Genotype", y="")+ theme(legend.position="none")
Tol_Sym_M4_SG.plot<-Tol_Sym_M4_SG.plot+labs(x="Site and Genotype", y="")+ theme(legend.position="none")
##Color scores
Tol_ScoreF_W2_SG.plot<-Tol_ScoreF_W2_SG.plot+labs(x="")+ theme(legend.position="none")
Tol_ScoreF_M1_SG.plot<-Tol_ScoreF_M1_SG.plot+labs(x="", y="")+ theme(legend.position="none")
Tol_ScoreF_M4_SG.plot<-Tol_ScoreF_M4_SG.plot+labs(x="", y="")+ theme(legend.position="none")
Tol_ScoreTP_W2_SG.plot<-Tol_ScoreTP_W2_SG.plot+labs(x="")+ theme(legend.position="none")
Tol_ScoreTP_M1_SG.plot<-Tol_ScoreTP_M1_SG.plot+labs(x="", y="")+ theme(legend.position="none")
Tol_ScoreTP_M4_SG.plot<-Tol_ScoreTP_M4_SG.plot+labs(x="", y="")+ theme(legend.position="none")
##Create Panel
Tolerance_fig<-plot_grid(Tol_Chl_W2_SG.plot, Tol_Chl_M1_SG.plot, Tol_Chl_M4_SG.plot,
Tol_ScoreF_W2_SG.plot, Tol_ScoreF_M1_SG.plot, Tol_ScoreF_M4_SG.plot,
Tol_ScoreTP_W2_SG.plot, Tol_ScoreTP_M1_SG.plot, Tol_ScoreTP_M4_SG.plot,
Tol_Sym_W2_SG.plot, Tol_Sym_M1_SG.plot, Tol_Sym_M4_SG.plot,
nrow=4, ncol=3, byrow=TRUE,
rel_widths=1,
rel_heights=c(1, 0.85, 0.85, 0.85),
labels=c("a", "b", "c",
"d", "e", "f",
"g", "h", "i",
"j", "k", "l"),
label_size=panel.lab.sz,
label_fontface = "bold")
##Save Figure
ggsave(filename="Figures/FigureS2_Tolerance_Across_Metrics.png", plot=Tolerance_fig, dpi=300, width=14, height=20, units="in")
##Save Figure
ggsave(filename="Figures/Figure4_Tolerance_Heatmap.png", plot=Tol_Pairs.plot, dpi=300, width=8, height=12, units="in")
##Combine Results Tables
TableS2_Tol.lm.res<-rbind(Tol_Chl_W2_lm.res, Tol_Chl_M1_lm.res, Tol_Chl_M4_lm.res,
Tol_Sym_W2_lm.res, Tol_Sym_M1_lm.res, Tol_Sym_M4_lm.res,
Tol_ScoreF_W2_lm.res, Tol_ScoreF_M1_lm.res, Tol_ScoreF_M4_lm.res,
Tol_ScoreTP_W2_lm.res, Tol_ScoreTP_M1_lm.res, Tol_ScoreTP_M4_lm.res)
##Add a Season Variable
TableS2_Tol.lm.res$Season<-ifelse(TableS2_Tol.lm.res$TimeP=="W2", "Summer1", ifelse(TableS2_Tol.lm.res$TimeP=="M1", "Summer2", ifelse(TableS2_Tol.lm.res$TimeP=="M4", "Winter" , NA)))
##Organize
names(TableS2_Tol.lm.res)
[1] "DF" "Sum.Sq" "Mean.Sq" "F.value" "p.value" "Predictor" "EtaSq"
[8] "Response" "TimeP" "Season"
TableS2_Tol.lm.res<-TableS2_Tol.lm.res[,c("TimeP", "Season", "Response", "Predictor", "DF", "Sum.Sq", "Mean.Sq", "F.value", "p.value", "EtaSq")]
#Round to 4 digits
TableS2_Tol.lm.res$Sum.Sq<-round(TableS2_Tol.lm.res$Sum.Sq, 4)
TableS2_Tol.lm.res$Mean.Sq<-round(TableS2_Tol.lm.res$Mean.Sq, 4)
TableS2_Tol.lm.res$F.value<-round(TableS2_Tol.lm.res$F.value, 4)
TableS2_Tol.lm.res$p.value<-round(TableS2_Tol.lm.res$p.value, 4)
TableS2_Tol.lm.res$EtaSq<-round(TableS2_Tol.lm.res$EtaSq, 4)
##Write Out Table
write.csv(TableS2_Tol.lm.res, "Tables/TableS2_Tolerance_LM_Results.csv", row.names=FALSE)
##Pairwise Results Table
TableS3_Tol.pairwise<-Tol_contrasts
##Add a Season Variable
TableS3_Tol.pairwise$Season<-ifelse(TableS3_Tol.pairwise$TimeP=="W2", "Summer1", ifelse(TableS3_Tol.pairwise$TimeP=="M1", "Summer2", ifelse(TableS3_Tol.pairwise$TimeP=="M4", "Winter" , NA)))
##Re-order by Seasonal Timepoint
TableS3_Tol.pairwise<- TableS3_Tol.pairwise %>% arrange(as.numeric(TimeP))
##Organize
names(TableS3_Tol.pairwise)
[1] "group1" "group2" "estimate" "SE" "df" "t.ratio" "p" "Sig"
[9] "Response" "TimeP" "Contrast" "Set" "Estimate" "Season"
TableS3_Tol.pairwise<-TableS3_Tol.pairwise[,c("TimeP", "Season", "Response", "Contrast", "Estimate", "SE", "df", "t.ratio", "p")]
TableS3_Tol.pairwise<-TableS3_Tol.pairwise %>% dplyr::rename( DF = df) %>% dplyr::rename( p.value = p)
#Round to 4 digits
TableS3_Tol.pairwise$Estimate<-round(TableS3_Tol.pairwise$Estimate, 4)
TableS3_Tol.pairwise$SE<-round(TableS3_Tol.pairwise$SE, 4)
TableS3_Tol.pairwise$t.ratio<-round(TableS3_Tol.pairwise$t.ratio, 4)
TableS3_Tol.pairwise$p.value<-round(TableS3_Tol.pairwise$p.value, 4)
##Write Out Table
write.csv(TableS3_Tol.pairwise, "Tables/TableS3_Tolerance_Pairwise_Results.csv", row.names=FALSE)